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PREFACE 


This publication consists of a series of project papers prepared by graduate stu- 
dents in computational fluid dynamics. The work was performed during the 1973-74 
academic year at Old Dominion University under the auspices of Professor Stanley G. 
Rubin, now with the Polytechnic Institute of New York. The class, taught at Langley 
Research Center, was comprised of Langley employees and contractors. Each paper 
briefly examines one of the numerical methods discussed during the lectures and 
applies the method to the Navier -Stokes equations governing incompressible flow in a 
driven cavity. Also, solutions obtained with a cubic spline procedure which was under 
investigation at the time are included. 

It should be emphasized that most of these studies were completed in a period of 
approximately 6 weeks by different researchers with backgrounds varying from fluid 
mechanics to applied mathematics. (All had considerable experience with numerical 
computations.) Therefore, this publication should not be used as a definitive comparison 
of the various numerical methods discussed, since this would require a more extensive 
and unified study. 

The results and discussion presented are instructive in that: 

(1) They demonstrate a number of numerical procedures applied by experienced 
numerical analysts who, however, were not expert in the use of these methods and were 
required to complete their analyses within a specified time period. 

(2) They include solutions for both primitive -variable and vorticity — stream - 
function systems of equations. 

(3) They present solutions for the driven cavity obtained with popular numerical 
procedures (such as alternating direction implicit, successive overrelaxation, and Mar- 
ker and Cell) as well as with several techniques that have not been applied extensively 
(e.g., direct Poisson solvers, spline methods, weighted differencing, artificial compres- 
sibility, and hopscotch). 

(4) They include some solutions for Reynolds numbers larger than those for which 
solutions have previously been reported in the literature. This publication is considered 
suitable for providing graduate students and practitioners of computational fluid dynam- 
ics with reference material covering a relatively broad spectrum of numerical proce- 
dures for solving the Navier -Stokes equations. 
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1. INCOMPRESSIBLE VISCOUS FLOW IN A DRIVEN CAVITY 


Stanley G. Rubin* and Julius E. Harris 
Langley Research Center 

SUMMARY 

In this paper the driven cavity problem is introduced, and the governing equations 
are given in terms of primitive and vorticity— stream -function variables. The results 
of the papers in this publication are discussed briefly. 

INTRODUCTION 

The flow of an incompressible viscous fluid in a cavity driven by a uniformly mov- 
ing boundary exhibits a number of complex fluid dynamic characteristics of interest 
including vortex development and boundary-layer development on the walls of the cavity 
(for high Reynolds numbers). Primarily because of the simple geometry, the incom- 
pressible driven cavity problem has received detailed attention in the literature as a test 
problem for evaluating numerical procedures for solving the Navier -Stokes equations. 
(See refs. 1 to 10.) Numerical results have been published for Reynolds numbers as 
large as 1000; however, these results indicate that most currently available numerical 
procedures deteriorate in accuracy and/or stability for Reynolds numbers much less 
than 1000. In most of these studies, numerical results were obtained for a Reynolds 
number of 100; therefore, the results for this Reynolds number should be used as a basis 
for the evaluation of the various numerical procedures presented in this publication. 

SYMBOLS 


h depth of cavity 

l width of cavity 

p pressure 

R Reynolds number, U l/v 

U velocity of moving upper surface 


♦Visiting professor at Old Dominion University, Norfolk, Va. Now professor of 
aerospace engineering, Polytechnic Institute of New York, Farmingdale, N.Y. 




u,v velocity component in x- and y-direction, respectively 

x,y rectangular coordinates defined in figure 1 

£ vorticity 

v kinematic viscosity 

p density 

if/ stream function 

Subscripts: 

t differentiation with respect to time 

x,y differentiation with respect to x or y 

STATEMENT OF PROBLEM 

In the present publication, incompressible viscous flow in a cavity driven by a uni- 
formly moving upper surface (see fig. 1) is predicted by several numerical techniques. 
Solutions are obtained for the Navier -Stokes and continuity equations written in either 
primitive or vorticity — stream -function variables. These equations are as follows: 

Primitive variables 

Continuity 


u x + Vy = 0 (la) 

x-momentum 

u t + uu x + vu y = -P x + R-^Uxx + Uyy) (lb) 

y-momentum 

V t + UV X + V Vy = -Py + R-^Vxx + Vyy) (l C ) 
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Vorticity— stream -function variables 


Vorticity 


^t + u ^x + vi =y ~ ^ (^xx + ?yy) 


(2a) 


Stream function 

^xx + ^yy = 5 (2b) 

All quantities are nondimensional. Velocities have been nondimensionalized with U, 
lengths with l, pressure with pU^, and time with 1/ U. 

The boundary conditions are presented in figure 1. In several papers contained 
herein the authors consider a plate moving from right to left so that u = i// y = -1 at the 
upper boundary; in addition, £ is defined as v x - u y . Direct comparisons with solu- 
tions for a plate moving with u = 1 are possible if u, v, and are replaced with 
-u, -v, and -\j/, respectively. 

DISCUSSION OF RESULTS 

In each paper included in the present publication the square driven cavity has been 
studied for a standard test case of R = 100. This particular case has received the most 
detailed attention in the literature. In addition to the standard test case, numerical 
results are presented in some papers for values of R of 10, 1000, and 5000. Some 
solutions are also given for rectangular cavities having width-depth ratios of 2 and 4. 
Solutions for R = 5000 were obtained with a central -difference alternating-direction 
implicit procedure for the vorticity equation and a successive overrelaxation or a direct 
Poisson solver for the stream function. These results for the driven cavity represent, 
to the best of the authors' knowledge, solutions for the highest Reynolds number currently 
available in the literature. Other results were obtained with hopscotch, Marker -and - 
Cell, spline integration, Crocco, weighted-differencing, Chorin, and Spalding procedures. 

The results presented in the papers indicate that convergence was more rapid with 
the vorticity— stream -function formulation. The accuracy of the primitive-variable solu- 
tions was extremely sensitive to the convergence tolerance assigned to the pressure sol- 
ver. As this criterion was refined, computational times escalated. 

Divergence form calculations were more accurate than nondivergence results, 
except for the nondivergence spline solutions. These were more accurate than any of the 
finite -difference results. This would reflect the higher order accuracy of the convection 
terms in the spline procedure. 
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It should be carefully noted that the computer codes developed for the procedures 
were not necessarily optimized; consequently, no comparisons among the papers are 
presented concerning computational efficiency (computer processing time and storage 
requirements). However, the numerical procedures are developed in sufficient detail so 
that they could be easily used by those interested in either making direct comparisons or 
applying the procedures to other problems. The intent of the authors is that the present 
publication should be an outline of the development of specific numerical procedures and 
their application to a standard test problem. From this viewpoint, the potential users of 
the various numerical techniques should be able to select the specific procedure which 
would best satisfy their particular computational requirements. 
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2. APPLICATION OF UPSTREAM WEIGHTED DIFFERENCING 
TO THE DRIVEN CAVITY PROBLEM 

Lillian R. Boney 
Langley Research Center 

SUMMARY 

The problem considered is steady plane flow of an incompressible viscous fluid in 
a two-dimensional square cavity in which the motion is driven by uniform translation of 
the upper wall which moves across the cavity from left to right at a constant speed of 1. 
Numerical solutions of the Navier -Stokes equations are obtained for a range of Reynolds 
number from 100 to 1000. Upstream weighted differencing is used when the local cell 
Reynolds number exceeds 2. 

INTRODUCTION 

It is well-known that the flow of viscous fluid can be treated successfully with the 
Navier-Stokes equations. Many authors have examined finite -difference solutions of the 
equations of steady-state motion and continuity in a two-dimensional domain. Some 
found that unless the mesh size was severely reduced, their computational procedures 
failed to converge when the Reynolds number became large. 

This paper treats steady plane flow of an incompressible viscous fluid in a two- 
dimensional square cavity in which the motion is driven by uniform translation of the 
upper wall which moves across the cavity from left to right at a constant speed of 1. 

(See paper no. 1 by Rubin and Harris.) The two momentum equations and the continuity 
equation, written in normalized conservative form in terms of the primitive variables 
(velocity components u and v and the pressure p), are used rather than the stream- 
function and vorticity transport equations. The pressure term is treated by a Poisson 
equation. The primitive-variable system involves two parabolic momentum transport 
equations and one elliptic equation with all Neumann boundaries. Upstream weighted 
differencing (ref. 1) for the velocity gradients is chosen for its good stability properties 
and the transportive property which finite -difference formulation of a flow equation pos- 
sesses if the effect of a perturbation in a transport property is advected only in the 
direction of the velocity. 
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SYMBOLS 


K constant in equation (5) 

p pressure 

p error in initial estimate of p 

p Q initial estimate of p 

R Reynolds number 

R c cell Reynolds number 

S = v2p 

t time 

At time step 

u axial velocity component 

v normal velocity component 

W weighting factor based on R c (see eqs. (7)) 

w factor in equation (22) 

x axial coordinate 

y normal coordinate 

Ax, Ay grid spacing in x- and y-directions 

fi = Ax/Ay 
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Superscripts: 


k index of iteration 

n index of time step 

Subscripts: 


i,j index of grid spacing in x- and y-direction, respectively 

ref reference 


GOVERNING EQUATIONS 


The governing equations are the two momentum equations and the continuity equa- 
tion. These are written as follows in forms of the primitive variables, velocity compo- 
nents u and v and pressure p, in a normalized conservative form: 


9u 9(u 2 ) + 3(vu) = _^P + 1 fd 2 u + 9 2 u \ 

9t 9x 3y 3x rU x 2 d y2J 

9v . 9(uv) 9(v 2 ) _ £P 1 /9 2 v 8 2 v\ 

3t 3x 9y 'ay Ry 9x 2 dy 2J 


9u 9v 
3x 3y 


= 0 


(3) 


where R is the dimensionless Reynolds number. The Poisson equation for pressure is 


V 2 p = S 


9u 9v 
3x " 9y 


(4) 


where 


P = P 0 - P 

and p 0 is the initial estimate of p. 
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FINITE -DIFFERENCE REPRESENTATIONS 


For numerical integration of the governing equations, the entire region of the cavity 
is divided into a finite number of grid points having coordinates x = i Ax and y = j Ay. 
The grid spacings in the x- and y-directions are denoted by Ax and Ay. An approxi- 
mate solution of the equations is obtained at these grid points at discrete times t n . The 
symbol t n denotes the time level after the nth time step. Practical stability imposes a 
restriction on the size of the time step At. Limited numerical experience suggests the 
following restriction: 


At = 


K 


2 

r 1 i i 

_L J 

u ul . 1 

R 

(Ax) 2 (Ay) 2 

T 

Ax 


M 

Ay 


(5) 


where K is generally taken to be 1. 

An explicit procedure was used to calculate quantities at time n + 1 in terms of 
their values at time n. The velocity gradients are calculated by an upstream weighted 
differencing method. (See refs. 1.) This method involves one-sided rather than space- 
centered differencing. Backward differences are used when the velocities are positive 
and forward differences are used when the velocities are negative. Thus the one-sided 
difference is always on the upstream side of the point at which the gradients of the veloc- 
ities are being evaluated. Upstream differencing is not stability limited by a cell Reyn- 
olds number 


R 


uy| + |vy| Ay 

VR 


( 6 ) 


A weighting factor W based on R c is introduced. For R c less than or equal to 2, 
the weighting factor becomes zero, and the differencing method becomes space -centered 
differencing. As R c approaches infinity, W approaches 1, and the differencing 
method becomes upstream differencing. That is, 


W = 0 


W = 1 


_ 2 _ 

R c 


(R c § 2) 

(Rc > 2) | 


(7) 
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Use of the weighting factor enables changes in flow direction to be accounted for, and 
adjustments are made with no discontinuity. For forward time differencing the equations 


are 






uft 1 - 

u^ 




9u - l >) 
at At 

hi 


(8) 


vPt 1 - 

vP- 




9V _ UJ 

9t At 



(9) 

The following terms are evaluated by central differencing: 



a 2 v _ v .i+l,j 

+ v i-l>j 

" 2v iJ 

(10) 


9x 2 

(Ax) 2 




9 2 v _ Vi »j +1 

+ v i,M 

- 2v i,i 

(ID 


ay 2 

(Ay) 2 




a 2 u _ Ui+1 »J 

+ u i-l»j 

- 2u i,j 

(12) 


9x 2 

(Ax) 2 



a 2 u _ U u+i 

+ u i,M 

- 2u i,j 

(13) 


8y 2 

(Ay) 2 



ap _ p i+l,j _p i-l,j 
ax 2 Ax 


(14) 


9p s p i,j+l ~ p i, j - 1 
9y 2 Ay 

The following terms are evaluated by upstream weighted differencing: 


(15) 


3u _ 1 
ax Ax 






W i+1J 






( 16 ) 
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I?" if + j^yj)( Ui,iVi,:l _Ui - 1 >j Vi - 1 ’j) + ?( 1 " |^j|J( u i+ 1 ,j v i+ 1 ,j - U u v u) 


+ ~^( u i + i,j v i+i,i - ^i-i,jVi_i,j) 


17 " i f (* + ]7j|]J( Vl *J ul »J ‘ V U -l u bj-l) + f (J - pT7|j( v i,i+l u iJ+l ' V U U U) 


W 4 v i, j 


+ ^T^( Vi ,j+ 1Ui ,j+ 1 “ v i,j-l u i,j-l) 


SPALDING METHOD FOR PRESSURE CALCULATION 


The Poisson equation for pressure (eq. (4)) with all Neumann boundaries on which 
values of the normal gradient 9p/9n are zero is solved by an iterative Liebman method. 
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For the initial iteration, p ; • is assumed to be zero at all grid points. For interior 
grid points, 


k+1 _ ~k 


PiJ = P 


w 


i,l 


s(i + £ 2 ) 


Pi+1,3 


+ Pui,1 


3 2 Py+i + - (Ax)2Si,j - 2(l + ^)pkj 


(22) 


where 



w = 1 


3u 9V 
ax ay 


An arbitrary limit is set to the number of iterations for the interior grid points, and the 
Neumann conditions are applied at the boundaries by setting the boundary points equal to 
the adjacent interior points. 

The pressure p re f is held at zero at a reference location chosen to be the center 
of the wall opposite the "driving" boundary. After an arbitrary number of iterations the 
pressure is corrected at all grid points by p = p Q - p^+l + p re f. The values of uP^ 

n 1 1 

and Vj j are corrected by setting 



The velocities u n and v n are replaced by u n+ ^ and v n+ *, and the procedure is 
repeated for the next time step. Numerical calculations were carried forward in time 
until a steady state was approached. Calculations were discontinued when the maximum 
values of ju n+ l - u n l and |v n+ l - v n | were less than 0.0005. Values of 


(u n +i _ u n )/u n+ l, (u n+1 - u n )/At, and ^ 


bl 


uPi 1 

- uP. 

.+ • 

v n+l 

- vP. 

1,1 

1,1 



bl 


were also com- 


puted for comparison. 


Approximately 0.25 second was required on the CDC 6600 computer for one time 
step and a 15 x 15 grid. 
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RESULTS 


Table I contains the number of time steps required to approach steady state at var- 
ious conditions. Reynolds numbers of 100, 400, and 1000 were chosen. Grid sizes of 
15 x 15, 29 x 29, and 57 x 57 and factors of 0.5, 1, 2, and 2.5 for the At restriction 
were used. Generally 10 iterations for p were calculated, and a check was made by 
using 5 and 20 iterations. 

Table II contains the numerical values of p, u, and v after 187 time steps for 
R = 100, a 15 x 15 grid, At = 0.0458, and 10 iterations of p. Figure 1 indicates the 
direction but not the magnitude of the velocity vectors and figure 2 is a contour plot of 
lines of constant pressure for the same case. Figures 3 and 4 show the velocity and the 
pressure profiles through the center of the vortex. 

The velocity vectors and velocity profiles through the center of vortex for R = 100 
compare well with those shown by Mills (ref. 2). The pressure distributions through the 
center of the vortex are similar but less smooth than those shown by Mills. The lines of 
constant pressure were compared with those shown by Kawaguti (ref. 3) for R = 64 and 
were less smooth. 

CONCLUDING REMARKS 

The method of upstream weighted differencing has been applied to the driven cavity 
problem. Results for the square cavity with a Reynolds number of 100 have indicated 
that the velocity information obtained is in good agreement with results obtained by other 
methods. The pressure information obtained does not agree as well and seems to be less 
reliable. 
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R 

100 

100 

100 

400 

1000 

1000 

1000 


Grid 

K 

(eq. (5)) 

At 

p iterations 

Time 

steps 

Convergence, 
u n+1 - u n 

15 x 15 

1.0 

0.0458 

10 

187 

0.48 x lO- 3 



.0458 

5 

220 

.45 x lO- 3 


1.0 

.0458 

20 

170 

.50 x lO' 3 

15 x 15 

.5 

.0229 

10 

220 

.49 x lO' 3 


2.0 

.0916 

10 

150 

.50 x lO" 3 


2.5 

.1144 

10 

20 

Nonconvergence 

29 x 29 

.5 

.0084 

10 

500 

.33 x 10-3 



.0168 

10 

360 

.50 x lO' 3 

15 x 15 


.0627 

10 

250 

.50x10-3 

15 x 15 

.5 

.0338 

10 

1730 

.46 x lO- 2 


1.0 

.0676 

10 

1170 

.87 x lO' 3 

29 x 29 

.5 

.0160 

10 

960 

.11 x lO' 2 


1.0 

.0321 

10 

600 

.13 x lO" 2 

57 x 57 

1.0 

.0146 

10 

240 

.20 x lO' 2 














TABLE II.- NUMERICAL VALUES OF PRIMATIVE VARIABLES AFTER 187 TIME STEPS 


[Each value computed at a grid point (i, j) ; columns correspond to i = 1, 2, 3, . . 15 

and lines to j 31 15, 14, 13, . . ., l] 








Pressure , 

p 







-.152 

-.152 

-.02 5 

-.049 

-.011 

-.022 

-.020 

-.024 

-.027 

- .021 

-.000 

.043 

. 101 

.264 

.264 

-.152 

-.152 

-.025 

-.049 

-.oil 

-.022 

-.020 

-.024 

-.027 

-.021 

-.000 

.043 

. 101 

.264 

.264 

-.064 

-.064 

-.036 

"*-.044 ' 

-.031 

-.036 

-.035 

-.037 

-.04 0 

- .034 

-.026 

.010 

.037 

. 127 

.12 7 

-.057 

-.05 7 

-.033 

-.039 

-.027 

-.031 

-.031 

-.034 

-.036 

-.031 

-.024 

.008 

. 033 

.095 

.095 

-.027 

-.02 7 

-.02 5 

-.028 

-.026 

-.029 

-.030 

-.033 

-.034 

-.030 

-.026 

-. 00 ? 

.009 

.041 

.04 1 

-.022 

-.02 2 

-.021 

-.021 

-.02 0 

-.022 

-.023 

-.024 

-.02 4 

-.020 

-.014 

.002 

.013 

.02 8 

. 028 

-.011 

-.OIL 

-.013 

-.014 

-.015 

-.016 

-.017 

-.017 

-.016 

-.012 

-. 008 

.003 

.010 

.014 

.014 

-.007 

-.00 7 

-.009 

-.009 

-.010 

-.010 

-.010 

-.010 

-.00 9 

-.005 

-.001 

.005 

.009 

. 008 

.008 

-.003 

-.003 

-.005 

-.005 

-.006 

-.006 

-. 006 

-.005 

-.004 

-.001 

. 002 

.005 

.007 

. 005 

.005 

-.001 

-.001 

-.004 

-.003 

-.004 

-.003 

-. 003 

-.002 

-.001 

.001 

.003 

.005 

.006 

.004 

.034 

-.000 

-.000 

-.002 

-.002 

-.002 

-.002 

-. 002 

-.001 

.000 

.002 

.004 

.004 

. 005 

. 003 

.003 

-.000 

-.00 0 

-. 00 ? 

-.001 

-.002 

-.001 

-.001 

.000 

.001 

.002 

.004 

.004 

. 005 

.00 3 

.00 3 

-.001 

-.001 

-.003 

-.002 

-.002 

-.002 

-.001 

-.000 

.001 

.002 

.003 

.004 

.004 

.003 

.003 

-.001 

-.00 1 

-.003 

-.003 

-.003 

-.002 

-.001 

0.000 

.002 

.003 

. 005 

.005 

.006 

.005 

.005 

-.001 

-.001 

-.003 

-.003 

-.003 

-.002 

-.001 

0 . 000 

.00 2 

.003 

.005 

.005 

. 006 

.005 

.005 

t 






Axial velocity 

» u 







1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1. 000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1 .000 

1.000 

0.000 

.130 

. 1 < J 3 

.2 89 

.340 

.402 

.446 

.484 

.512 

.526 

.517 

.487 

. 382 

.296 

0.000 

0.000 

-.040 

-.020 

.023 

.060 

.112 

. 150 

. 190 

.216 

.236 

.225 

.206 

.129 

.084 

0.000 

0.000 

-.056 

-.065 

-.053 

-.041 

-.012 


.032 

.053 

.067 

.030 

.011 

-.029 

-.921 

0.00 0 

0.000 

-.030 

-.048 

-.061 

-.068 

-.065 

-.063 

-.060 

-.062 

-.067 

-.081 

-.086 

-. 091 

-.059 

0.000 

U .000 

-.02 3 

-.044 

-.065 

-.081 

-.092 

-. 103 

-.111 

-.122 

-.131 

-. 140 

-.133 

-.112 

-.065 

0.000 

0.000 

-.016 
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Figure 1.- Direction of velocity vectors. R = 100; 15 x 15 grid. 
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Figure 2.- Contour plot of lines of constant pressure. R = 100; 15 x 15 grid. 
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X 

(a) Horizontal pressure profile. 

Figure 4.- Pressure profiles through center of vortex. R = 100; 15 x 15 grid. 
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3. SOLUTION OF THE INCOMPRESSIBLE DRIVEN CAVITY PROBLEM 
USING THE CROCCO METHOD 

Jerry N. Hefner 
Langley Research Center 

SUMMARY 

The Crocco method has been used to solve the steady Navier -Stokes equations for 
two-dimensional flow in a square driven cavity for a Reynolds number of 100. The 
The Spalding SIMPLE procedure was used to solve the pressure -linked equations. 
Results for velocity and pressure were found to be inaccurate, although the overall 
recirculation and vortex center were correctly predicted. 

INTRODUCTION 

During the past decade, the evolution of large high-speed computers has allowed 
the development of numerous explicit and implicit finite -difference schemes for solving 
the steady Navier -Stokes equations. In general, these numerical procedures have not 
been rigorously tested on a variety of problems. Consequently, the value of a given 
method for solving a multiplicity of fluid dynamics problems often has not been ascer- 
tained. In this paper, results of a numerical study to predict the velocity profiles and 
pressure distribution in the simple two-dimensional driven cavity with width -depth ratio 
of 1 (square cavity) and Reynolds number of 100 are presented. (See paper no. 1 by 
Rubin and Harris for details of driven cavity problem.) The method used is that attri- 
buted to Crocco (ref. 1) with the Spalding SIMPLE procedure (ref. 2) for solving the 
pressure -linked equations. 


SYMBOLS 


l width of cavity 

A l spatial increment in x and y 

p pressure 

p error in pressure (eq. (9)) 

p 0 initial estimate of p 
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q 


actual velocity vector 


q* calculated velocity vector (eq. (9)) 

R Reynolds number, ^ 

v 

t time 

At time step 

U velocity in x-direction at edge of cavity 

u velocity in x-direction 

v velocity in y-direction 

x,y axial and normal coordinate of cavity, respectively 

T variable used in equations (5) to (8) 

v viscosity 

p density 

a) variable in successive overrelaxation method given in equation (12) 

Superscripts: 

k iteration number 

n time step number 

Subscripts: 

i step in x-direction 

j step in y-direction 
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t differentiation with respect to time 

x,y differentiation with respect to x or y 

Bars over symbols denote dimensional quantities. 

METHOD OF ANALYSIS 

The nondimensional viscous incompressible flow equations in primitive variables 

are 


U X + Vy = 0 


(1) 


u t + uu x + vu y + P x = l( Uxx + u yy ) 

V t + UV X + Wy + P y = i( VxX + Vyy) 

The variables have been nondimensionalized with respect to U and F; that is, 



U 




l 



( 2 ) 
: ( 3 ) 


(4) 


The boundary conditions are as follows: 
at y = 0 u = v = 0 

at y = 1 u = 1 and v = 0 

at x = 0 u = v= 0 

at x = l u = v = 0 

The Crocco method (ref. 1) was applied to the x- and y -momentum equations (eq. (2) 
and (3), respectively) to solve for velocities u and v. The method presented by Crocco 
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is actually an iteration procedure which avoids the difficulties associated with numerical 
solution of the steady Navier -Stokes equations by considering the steady flow as the 
asymptotic form of a time -dependent flow. These transient flow calculations, although 
both inaccurate and inconsistent, are used only as an iteration procedure and permit 
freedom in choosing numerical differencing schemes which best fit the desired stability 
requirements. In the Crocco method, the convective terms are redefined in terms of 
their values at time steps n and n - 1 by means of a variable T which can be 
chosen to account for various differencing schemes (e.g., for r a 0, the differencing 
scheme is simply forward time and centered space, whereas for T = 1/2, the method 
becomes the well-known Adams -Bashforth method). As an illustration of the redefinition, 
the convective terms in the x -momentum equation (2) become 

uu x = (1 + T)u n u x n - ru n-1 u§ _1 (5) 

vuy = (1 + r)v n Uy n - rvH-iup- 1 (6) 

Modified forward time and centered space differencing (i.e., the values of uy and 
vy in the spatial derivatives are evaluated at time n + 1 instead of n) is applied to 
the appropriate derivatives. Values of r between 0 and 1.0 are chosen depending on 
the desired stability requirements. (For stability, the Crocco method requires that 
T > 1/2.) 

Application of the Crocco method to equations (2) and (3) gives the following differ- 
ence equations: 

x -momentum 


1 + 


4 At 

r(mV 


u: 


n+1 At 


i,J 


R(AZ) 


iHki ♦ j + "y+i - <m) - - <u) 
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y-momentum 



4 At 
R(Af) 2 


,n+l 

i, 


At 

R(AZ) 2 



+ v 


U+l 


Km) 


(i + r) u fj it 


2 M 


(< 


ri+ij 



(i + nvfj At - 

l »J -M 1 - i -yP- 


2 M 


rdM At 
l ?J 

2 A£ 


- *-£i 


r'-y 1 At 

2 M 


H*»j+1 " V U- J ) " FazRj+I " p i,j-l) 


(8) 


with the spatial increments in x and y being equal and defined by A l. 

The Spalding SIMPLE procedure (ref. 2) was used to solve the pressure -linked 
Poisson equation. The SIMPLE procedure is: 

(1) Assume an initial value p Q for pressure. 

(2) Calculate u and v based on the assumed p Q by any method. 

(3) Let the calculated velocity vector be q*. 

(4) Assume that the actual velocity vector is 


q = q* + Vp 


(9) 


The divergence of both sides of equation (9), along with the continuity equation, gives 


V 


q = 0 = V • q* + V 2 p 


(10) 


The surface normal pressure gradients are assumed to be zero. 

(5) The actual pressure p is then defined as 

P = P 0 + P (11) 
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The successive overrelaxation method is applied to equation (10) and the differ- 
encing gives 


pk+l,n+l = pk,n+l + w pk,n+l + pk+l,n+l + pk,.n+l -k+i n+i 


4 'k,n+l Ai/( n+1 

4p i,j ' TH+l.J 


n+1 n+1 

u i-l,j + v i,j+l 


v 


n+1 \ 

y-ij 


RESULTS 


( 12 ) 


The results presented in table I and figures 1 and 2 are representative of the more 
accurate solutions obtained in the present study for the square cavity flow problem. The 
solution presented was obtained for At = 0.035, A l = 0.0714, optimum = 1.5456, and 
T = 0.6. Comparing the velocity profiles at the vortex center with those of reference 3 
shows the inaccuracy of the present solution. The present predictions for velocities u 
and v disagree with reference 3 by as much as 10 to 12 percent and 30 to 35 percent, 
respectively. Furthermore, examination of table I indicates significant oscillations in 
the calculated pressures throughout the cavity. Despite these inadequacies the present 
solution does predict the correct overall recirculation and vortex center in the cavity. 

(See fig. 1.) 

Attempts to improve the accuracy and convergence of the steady-state solution 
were unsuccessful in computing times up to 3000 sec on the CDC 6600 computer. Vary- 
ing the time step by a factor of 2, varying r from 0 to 1.0, and using forward time and 
centered space differencing instead of modified forward time and centered space differ- 
encing on the viscous terms did not significantly affect the accuracy or convergence of the 
present solutions. In reference 3, Mills obtained accurate solutions by solving the equa- 
tions for the cavity flow in terms of stream function and vorticity; therefore, the need to 
calculate pressure directly was eliminated. The present solutions using primitive vari- 
ables require the pressures to be calculated, and hence the accuracy of the overall solu- 
tion is dependent on the accuracy of the pressure calculation. Therefore, the relatively 
poor accuracy of the present solutions and the oscillations in predicted pressures are 
probably attributable to use of primitive variables and inaccurate solutions to the 
pressure -linked equations. 


CONCLUDING REMARKS 

Solutions of the steady Navier -Stokes equations using the Crocco method along witu 
the Spalding SIMPLE procedure have been obtained for two-dimensional flow in a square 
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driven cavity. The results for the axial and normal velocities did not agree with those 
of another accurate method, and significant pressure oscillations were noted. These 
inadequacies were probably due to use of primitive variables and to inaccurate solutions 
to the pressure-linked equations. The solution did correctly predict the overall recir- 
culation and vortex center. 
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TABLE I.- THEORETICAL PRESSURE DISTRIBUTION IN SQUARE CAVITY 


Pressures at x of - 



0 

0.071 

0.143 

0.214 

0.286 

0.357 

0.429 

0.500 

0.571 

0.643 

0.714 

0.786 

0.857 

0.929 

1.000 

1.000 

0.581 

0.570 

0.622 

0.626 

0.638 

0.640 

0.628 

0.637 

0.613 

0.648 

0.631 

0.731 

0.758 

0.972 

1.000 

.929 

.562 

.520 

.628 

.604 

.624 

.628 

.598 

.636 

.565 

.655 

.554 

.749 

.645 

1.058 

.998 

.857 

.612 

.619 

.626 

.627 

.624 

.627 

.610 

.591 

.619 

.590 

.661 

.668 

.668 

.807 

.863 

.786 

.614 

.584 

.623 

.606 

.616 

,619 

.595 

.623 

.570 

.633 
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.714 
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.611 
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.625 
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.659 

.652 
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.700 

.571 

.653 

.644 

.647 

.643 

.642 

.640 

.635 

.635 

.631 

.634 

.639 

.641 

.661 

.639 

.670 

.500 

.650 

.625 

.644 

.628 

.638 

.632 

.632 

..637 

.630 

.644 

.639 

.653 

.656 

.549 

.666 

.429 

.661 

.651 

.656 

.650 

.653 

.649 

.651 

.648 

.651 

.648 

.656 

.648 

.661 

.637 

.659 
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Figure 1.- Theoretical flow -field predictions for square cavity. 
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4. SOLUTION OF THE DRIVEN CAVITY PROBLEM IN 
PRIMITIVE VARIABLES USING THE SM AC METHOD 

Richard S. Hirsh 
Langley Research Center 

SUMMARY 

The equations of two-dimensional incompressible fluid flow in the driven cavity 
have been solved in terms of primitive variables. The Simplified Marker -and-Cell 
method was used with forward time and centered space (FTCS) differencing and with 
modified forward time and centered space (MFTCS) differencing. The problem could 
not be solved for a Reynolds number of 100 with FTCS differencing because of stability 
limitations. For a Reynolds number of 10, the results were found to be inaccurate 
because of the inability of the pressure solver to converge to small tolerances. The 
MFTCS differencing seemed to relax the stability restriction on time step, but the limit 
on cell Reynolds number was retained. A lower order boundary condition on pressure 
resulted in marginal improvement. 


INTRODUCTION 

The equations of two-dimensional incompressible fluid flow can be integrated by 
either of two general procedures. The equations can be utilized as they are, written 
expressing conservation of mass and momentum in terms of axial and normal velocity 
and pressure, that is, the so-called primitive variables (PV); or the continuity equation 
can be satisfied by a stream function, and the pressure can be eliminated by taking the 
curl of the momentum equations to yield the vorticity— stream -function (VSF) system. 

The VSF system is attractive since it reduces the number of equations from three 
to two and eliminates the troublesome continuity equation. This gain is offset somewhat 
by the difficulties encountered when dealing with the unknown boundary conditions on 
vorticity. In some problems, such as aerodynamic flows, for which knowledge of the 
pressure is desired, the VSF system must be supplemented by additional relations for 
the pressure. The PV system, on the other hand, yields the pressure directly and has 
little difficulty with the boundary conditions. The present study considers the solution 
of the driven cavity problem with a method which solves the PV system, and results are 
compared with a solution of the VSF system. See paper no. 1 by Rubin and Harris for 
details of the driven cavity problem. 
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SYMBOLS 


n unit surface normal 

p pressurelike variable 

q correct velocity vector 

q* approximate velocity vector 

R Reynolds number 

R c cell Reynolds number 

t time 

At time step 

u axial velocity 

v - normal velocity 

x,y axial and normal coordinate, respectively 

Ax, Ay spatial increments in x- and y-directions 

A dilatation 

e convergence factor in SOR pressure solver 

Superscripts: 
n time level 

* intermediate time level 
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Subscripts: 


i,j index of grid points in x- and y-direction, respectively 

x,y,t differentiation with respect to x, y, or t 

Abbreviations: 

ADI alternating direction implicit 

FTCS forward time and centered space 

MAC Marker and Cell 

MFTCS modified forward time and centered space 
SMAC Simplified Marker and Cell 

SOR successive overrelaxation 

METHOD OF SOLUTION 

The momentum equations in two-dimensional incompressible flow, written in 
dimensionless divergence form, are (ref. 1) 

Uf + (u 2 ) x + (uv) y = -P x + ^(Uxx + Uyy) (1) 

V t + (uv) x + (v 2 ) y = -p y + i(vxx + Vyy) r (2) 

The continuity equation is 

A = u x + Vy = 0 (3) 

These three equations constitute the governing equations of the driven cavity problem. 

The method of solution used is the Simplified Marker -and-Cell (SMAC) method 
originally proposed by Amsden and Harlow. (See ref. 2.) Although this method does not 
yield the true pressure from the solution, it gives an indication of the accuracy which can 
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be achieved with the equations in terms of primitive variables since the velocity is cor- 
rect, as pointed out in reference 2. 

The methodology of the SMAC method is as follows: 

(1) Integrate the momentum equations one time step to generate an approximate 
vector solution q*. 

(2) Knowing this to be incorrect since in general, it does not satisfy zero dilatation, 
assume that a pressurelike variable p corrects q* to the actual value q such that 
div q = 0. Thus, say 

q = q* + grad p (4) 

(3) Since div q = 0, equation (4) gives 

v2p = _div q* (5) 

(4) Since q = q* on the boundaries, equation (4) gives for boundary conditions, 


= o 

an 


( 6 ) 


(5) Once the Poisson equation (5) for the "pressure " is solved, the correct value 
of q at this particular time step can be found from equation (4). Then, the whole pro- 
cess is begun again at step (1) and continued in time until two consecutive time steps 
have velocities which differ by no more than some small convergence factor. 

Note that the variable p is not the real fluid pressure, since the boundary condi- 
tion (eq. (6)) is not correct. The true normal conditions on the pressure can be derived 
from equations (1) and (2) written at the wall. 

The spatial differencing used in the numerical solution of equations (1), (2), and (5) 
is based on the staggered grid originally used in the MAC method (ref. 1), that is, 
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With forward time and centered space (FTCS) differencing, equations (1) and (2) are 
written 


u f+l/2,j ~ u f+l/2,j . Wi+jU ~ ( u \j [ (uv) f + l/2,j+l/2 - ^ uv ^+l/2,j -1/2 _ 

Ay 


At 


Ax 


pp . . - pp . 

1+1, J 1,] , 1 

u ?+3/2,j - 1 

L 2u F+l/2,jj 

f + u ?-l/2,j 

Ax R 


(Ax)^ 



4 


u i+l/2,j+l - { 2u i+l/2,j} +u i+l/2,i-l 


(Ay)' 


(7) 


'i*,j+l/2 - v i!j+l/2 ( uv >r+l/2,j+l/2 ~ «-l/2, j+ l/2 , ^U+l " = 

Ay 


At 


Ax 



where 


U i+l,j " |( u i+3/2,j + u i+l/2,j) 


u i,j = |( U i+V2,j +u i-l/2,j) 


37 



V i>j + 1 2( V U +3 / 2 + V i>j + l/2) 


( uv )i+l/2,j+l/2 = J(ui+l/2,j+l + u i+ l/2,j) ( v i+ l,j +1/2 + v i,j+l/2) 


( uv )i+l/2,j -1/2 = |( u i+l/2,j + u i+l/2,j -l) ( v i+l,j -1/2 + v i,j-l/2) 
( uv )i-l/2,j+l/2 = |( u i-l/2,j+l + u i-l/2,j)( v i,j+l/2 + v i-l,j+l/2) 


The difference equation used in the solution of the Poisson equation is 


p i+l,j ~ 3Pj,j + Pi-1, j p i,j+l _ 2p i,j + p i,j-l _ 


(Ax) 4 


(Ay) 4 


n * it* it* 

U i + l/2,j ~ *4-1/2,] V i,j+l/2 - v i,j -1/2 
Ax Ay 


(9) 


A simple successive overrelaxation (SOR) routine using optimal relaxation parameters 
was written to solve this equation rather than an alternating-direction implicit (ADI) tech 
nique since the efficiencies of these two methods are equivalent if no acceleration is used 
with the ADI method. The SOR routine was assumed to be converged when two succes- 
sive iterations differed by no more than a prescribed convergence factor. 

The use of FTCS differencing imposes stability restrictions on the time step, as 
well as on the cell Reynolds number. These three restrictions, given in reference 1, 
are (with equal mesh spacing) 


/3 = — < I 

R(Ax) 2 4 


Cx + C y= (|U| + |v|)|is l)> 


( 10 ) 


Rc = (M + |v|)R Ax = 4 
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If the maximum cell Reynolds number is assumed to occur at the upper moving wall of 
the driven cavity, the maximum flow Reynolds number which could be expected to be 
stable for a 15 x 15 grid is 56. Thus the test case of R = 100 could not be calculated 
with FTCS differencing. A value of R = 10 was chosen as the test case for runs using 
FTCS differencing in order to be on the conservative side of the stability limit. At these 
low values of R, the /3 -restriction in equation (10) is more severe than the c -restriction, 
so the test case was first run at an again conservative value of /3 = 0.1. After these 
runs proved to be successful, the values of /3, R, Ax, and the convergence factor in the 
SOR solver were varied. The results of these runs are described in the next section. 

In an attempt to overcome the restriction on /3, a modified time integration scheme 
(MFTCS differencing) was used. In this approach the braced terms of the viscous 
expressions in equations (7) and (8) are evaluated at time * instead of time n. This 
differencing scheme has been analyzed by Rubin and Lin (ref. 3) in their appendix for the 
one -dimensional case. Generalizing to two dimensions gives for stability when R c > 4, 

i (11) 

Rc - 16 

and gives unconditional stability for /3 when R c ^ 4. Thus the restrictions (eq. (10)) 
are somewhat relaxed. Runs duplicating the calculations with FTCS differencing were 
made with MFTCS differencing and are discussed in the following section. 

DISCUSSION OF RESULTS 

The standard problem with R = 10 was first run on a 15 x 15 grid with a crude 
convergence factor in the SOR pressure solver (e = 10 - 3). The flow in this and all other 
runs was assumed to be converged to a steady state when the velocities at two consecu- 
tive time steps differed by less than 10~6. As with all the results to follow, this first 
run gave a flow pattern which was qualitatively correct, a one -cell circulation within the 
cavity. (See fig. 1.) However, the final answer generated values for the dilatation 
A = div q which were not close to zero. Values in excess of 0.2 were generated close 
to the boundaries, and the entire A field did not give the appearance of a quantity which 
should be identically zero. 

In order to have a standard of comparison, the ADI solution of the next paper by 
Morris (paper no. 5), which agreed exactly with the published results of Mills (ref. 4) 
for R = 100 (and was therefore presumed to generate the correct result) was run for 
R=10 on a 15 x 15 grid. The SMAC results and the ADI results for the axial velocity 
u at the first row of grid points from the moving wall are shown in figure 2. The differ- 
ence between the maximum values of u is about 7 percent. It was felt that the crude 


39 



convergence for the pressure might lead to inaccuracies in the velocities (through the 
pressure gradients in the momentum equations), so the convergence factor for pressure 
was lowered an order of magnitude to e = 10-4, and the identical problem was run again. 
The computation time was increased by a factor of 2. 5, but the converged dilatation field, 
listed in table I, was considerably improved. There was a favorable change in all the 
velocities as shown in figure 2. Now the difference between maxima was only 
4.5 percent. 

Since the accuracy of the pressure calculation seemed crucial, the SOR conver- 
gence factor was decreased another order of magnitude to e = 10'^. With this value the 
SOR pressure solver would not converge. It was then decided to use e = 10-5 but to cut 
off the iteration at a large arbitrary limit (1000 iterations) if the convergence was not 
fulfilled and proceed to the next time step. The result was a significantly more uniform 
dilatation field. However the typical value of dilatation computed, 0.0040, was no better 
than the value for e = 10 -4 (table I). Thus, it seems that in addition to being the most 
time-consuming segment of the calculation, the method used for the pressure solver con- 
trols the overall accuracy of the computation. A plot of the velocity vector field for 
e = 10 “4 is presented in figure 1. The calculated velocities through the vortex center 
are compared with those of paper no. 5 in figure 3. 

Attempts were made to solve the standard problem (i.e., R = 100) with FTCS dif- 
ferencing. As expected, the calculation diverged. At R = 50, the solution also diverged, 
but at R = 40, the computation did converge and gave a dilatation field comparable to 
that computed for R = 10. This behavior seemed to confirm the stability criteria of 
equation (10). 

The MFTCS differencing was tested on the 15 x 15 grid: For R = 10, time steps 
20 times larger than those used with FTCS differencing and 8 times larger than those 
indicated by the stability restriction (eq. (10)) were used. The results with MFTCS and 
FTCS differencing were identical. This was also true for R = 40. Although better sta- 
bility was expected with MFTCS differencing (see eq. (11)), the solution for R = 50 
diverged, even for a very small /3, No explanation for this behavior could be found. 

The MFTCS method did require less total computer time than the FTCS computations, so 
all further runs were made with MFTCS differencing and a time step larger than the |3 
stability limit of equation (10). 

In an attempt to improve the accuracy for R = 10, the grid size was halved, and 
the computation was performed on a 29 x 29 grid. With e = 10 _ 3, the results (at com- 
parable points) were essentially the same as those on the 15x15 grid at e = 10-4. 
Computer time required was 6 times as much. Decreasing the convergence to e = 10 -4 
for the fine grid did affect the velocity field at a location y = 1/14 away from the upper 
wall as shown in figure 2. Now the difference between the maxima of the ADI results and 
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the SMAC results was only 2.3 percent; however, the dilatation remained around 0.1 for 
most points. It must be remembered, however, that these comparative points in the 
29 x 29 grid are no longer the ones adjacent to the boundary walls, so any wall effects 
felt by the first line of the 15 x 15 grid may be damped in the 29 x 29 grid by intervening 
points. 

Some simple experiments with the normal gradients used at the boundaries in the 
pressure solver were performed to see the. effects of different orders of accuracy. The 
previous results were all obtained with a second-order accurate end difference for the 
gradient; that is, 

9p _ 3 Pi,l - 4 Pi,2 + Pi, 3 _ Q 

3y i( l 2 Ay 

A first-order accurate difference at the wall is 

8p _ Pj,2 -Pl,l 

8y 1,1 iy 

which was inserted into the SOR routine, and all the previous runs were repeated. No 
substantial changes due to this new boundary condition were noticed. The velocities 
seemed higher in some instances, but in others the dilatation field seemed a bit lower. 

CONCLUDING REMARKS 

The Simplified Marker -and-Cell method using primitive variables did not give 
accurate results for flow in the driven cavity because of the inability of the pressure sol- 
ver to converge to small tolerances. Perhaps, instead of correcting the velocity vector 
only once at each time step, additional pressure corrections could be made until the dila- 
tation is reduced to a small level. Then a new time step could be calculated. This might 
also produce a more consistent scheme. 

Modified forward time and centered space differencing seems to relax the stability 
restriction on the time step but still retains a limit on cell Reynolds number. The lower 
order boundary condition on the pressure gave marginal improvement in the dilatation, 
but the computed values of velocity were slightly higher. 
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TABLE I. - DILATATION AT INTERIOR MESH POINTS FOR R = 10, e = 10-4, AND 15 X 15 GRID 
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5. SOLUTION OF THE INCOMPRESSIBLE DRIVEN CAVITY PROBLEM 
BY THE ALTERNATING -DIRECTION IMPLICIT METHOD 

Dana J. Morris 
Langley Research Center 

SUMMARY 

A second-order accurate alternating-direction implicit (ADI) method has been used 
to solve the two-dimensional incompressible Navier -Stokes equations for flow in a square 
driven cavity. Calculations were made at a Reynolds number of 100 with equally spaced 
15 X 15, 17 x 17, 33 x 33, and 57 x 57 grids and at a Reynolds number of 10 with a 
15 x 15 grid. The ADI results agreed well with other solutions. Choosing a criterion 
for convergence which was too strict was found to result in no steady-state solution being 
reached. A study of the maximum allowable time step for a Reynolds number of 100 
indicated that time steps many times larger than those used for most consistent explicit 
methods could be used for the ADI method. 

INTRODUCTION 

A second-order accurate efficient implicit method is used to solve the two- 
dimensional incompressible Navier -Stokes equations for flow in a square driven cavity. 
(See paper no. 1 by Rubin and Harris for details of the driven cavity problem.) This 
problem was previously solved by Mills (ref. 1) using an explicit method. 

Alternating -direction implicit (ADI) methods were introduced by Peaceman and 
Rachford (ref. 2) and take advantage of a splitting of the time step for multidimensional 
problems to obtain an implicit method which requires only the inversion of a tridiagonal 
matrix. The expected unconditional stability of the method allowed larger time steps 
than explicit methods without the complexity of a full matrix reduction routine required 
by one -step implicit methods. 




SYMBOLS 

K 

constant 


M 

>> 

< 

ii 

X 

< 

ii 


N 

number of grid points in 

one direction 
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R 


Reynolds number 


At time step 

u axial velocity component 

v normal velocity component 

x,y axial and normal coordinate, respectively 

Ax, Ay spatial increments in x- and y-directions 

£ vorticity 

u kinematic viscosity 

if/ stream function 

Superscripts: 

n index of time step 

* intermediate time level 

Subscripts: 

i,j index of grid point in x- and y-direction, respectively 

t differentiation with respect to time 

w value at wall 

x,y differentiation with respect to x or y 

METHOD OF SOLUTION 

The nondimensional equations for viscous incompressible two-dimensional flow in 
vorticity — stream -function form which describe the flow field under investigation are 

?t = _u ^x ” v ^y + ^xx + v ^yy 
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(1) 



V 2 l p= C 


( 2 ) 


and u = \py, v = -xp x , and v = 1/R. 

Applied to the vorticity transport equation (eq. (1)), the Peaceman-Rachford (ref. 2) 
alternating -direction implicit (ADI) method advances one time level in the following two 
steps. 

Step 1 


iiiLliid = _ u? * 

At/2 


V?y n + 


+ vK 


n 

yy 


Step 2 


v.n+1 _ 


= -u - v?J +1 


+ K 


n+1 

yy 


The value has no physical meaning. 

The method applied to the linear equation has a formal error of order 

o o o 

O (At) ,(Ax) ,(Ay) . It is unconditionally stable and consistent. The full second-order 
accuracy of the method can be deteriorated by the nonlinear terms which should be eval- 
uated as u* and v n for step 1 and as u* and v n+ l for step 2. If the old values 

[ 2 2 

At, (Ax) , (Ay) . Briley (ref. 3) calculated 

u n and v n and, from the previous values u 11 "* and v n_ *, linearly extrapolated for- 
ward to u* and v*. This procedure is stable and appeared second-order accurate; 
however, additional storage is required for Another procedure (refs. 4 and 5) is 

to iterate on the entire time step, either to convergence or for a predictor-corrector 
method. In either case the error is o[(At) 2 ] and additional storage is required for 
i// n+ *. It is also possible to calculate i p* after step 1 and to use u* and v* in 
step 2. Aziz and Heliums (ref. 5) examined these alternatives and found iteration on the 
full two steps to be the most accurate. The iteration required to obtain full second-order 
accuracy of the nonlinear terms may not be undesirable since some iteration is an 
advantage in achieving accurate boundary values for £. 

Although a linear Von Neumann analysis shows the ADI procedure to be uncondition- 
ally stable, a survey of the literature indicates that the procedure may not be uncondition- 
ally stable for flows with high Reynolds numbers because of the At time lag of £ w . (See 
ref. 6.) The degree of convergence for £ w to obtain stability is problem dependent. 

For large At, convergence may be prohibited by nonlinear effects. Also, the program - 
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ing involved with complex boundaries effectively restricts the use of this method to rec- 
tangular regions. 

The tridiagonal nature of the matrix to be solved, although much simpler than that 
in fully implicit one -step methods, requires diagonal dominance if error buildup is to be 
avoided in the matrix inversion. This is equivalent to imposing a restriction on cell 
Reynolds number. 

In the present study, the ADI procedure was applied to the vorticity transport equa- 
tion. Central differencing was used with equal mesh spacing and nonlinear terms were 
lagged a full time step. After completing a full time step calculation on vorticity, the 
stream-function equation was solved by a successive overrelaxation (SOR) routine and 
the boundary conditions were updated. (The SOR routine was identical with the one used 
in the previous paper by Hirsh (paper no. 4).) The vorticity was then calculated for a 
new time step, rather than being iterated at each time level to achieve a solution which 
was second -order accurate in time, since the interest in this study was in the steady- 
state solution. 


BOUNDARY CONDITIONS 


The no-slip condition requires that both u and v (the normal and tangential 
gradients of are zero at the stationary walls whereas at the moving wall u is unity. 
These conditions are treated by taking a row of image points at a distance AZ = Ax = Ay 
outside the boundaries. These values can be related to the interior points by taking 
derivatives at the boundaries. Thus, the boundary conditions for £ are: 

on the stationary walls, 

t° = ^interior 


at the moving wall, 

_ 2(AZ + ’/'interior) 

(Ai)2 

The boundary value of i/' is zero on all boundaries. Initial conditions are zero every- 
where except on the moving wall. 
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RESULTS 


These calculations were made at a Reynolds number of 100 for a square N x N 
cavity for N = 15, 17, 33, and 57 and at a Reynolds number of 10 for N = 15. For 
N = 17, the problem was also computed with the conservation form of the equations and 
the values were found to be much more accurate, although the computer time to achieve 
a converged solution was approximately 2 times longer. As shown in table I, the maxi- 
mum value of the stream function increases with the number of grid points toward the 
value 0.101 published by Burggraf (ref. 7). The vorticity £ at the midpoint of the mov- 
ing wall is also shown. For N = 15, the solution was identical with the one published by 
Mills (ref. 1) with a maximum stream function of 0.08742. Machine plots of the velocity 
vectors and the streamlines of this flow field are shown in figure 1. The calculated 
values of vorticity £ and stream function \J/ for N = 57 are presented in table II. 

In figure 2 the axial component of velocity u is depicted along a vertical line passing 
through the vortex center, and the calculated values are shown in table m for the three 
cases, N = 15 and 57 in nondivergence form and N = 17 in divergence form. Attempts 
to increase the accuracy of the nonconservation form by increasing the number of nodal 
points resulted in solutions which did not converge because of a convergence criterion 
which was too strict. The percentage difference between time levels n + 1 and n was 
used to check convergence. For N = 15 and 17, solutions were achieved with the con- 
vergence requirement set at 0.00005. For N = 33, this requirement was repeatedly 
made less stringent and was set at 0.001 before a converged solution was reached. Solu- 
tions for N = 57 did not require further weakening of this convergence criterion. This 
study was made for a time step At = Ax, the Courant-Friedrichs-Lewy condition for 
explicit methods. Decreasing the Reynolds number to 10 for N = 15 also forced the 
convergence requirement to be set at 0.001 and the time step to be decreased. No 
attempt was made to find the maximum allowable time step for this Reynolds number. 

For a Reynolds number of 100, however, such a study was made with the nonconservation 
form of the equations and the results are shown in figure 3. For N = 15, no solution 
would converge for At greater than 7 Ax. An N of 17 required that At be lowered 
to 5 Ax. Increasing N to 33 with At = 5 Ax resulted in the solution iterating in the 
SOR routine until the time limit on the computer was reached. Lowering the time step to 
3 Ax allowed the solution to converge. For N = 57, convergence in the SOR routine 
required that At = Ax. Thus the maximum allowable time step appears to vary not only 
as a function of Ax but also with the number of grid points used. 
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CONCLUSIONS 


An alternating -direction implicit (ADI) method has been used to solve the two- 
dimensional incompressible Navier -Stokes equations for flow in a square driven cavity. 
The following conclusions may be drawn: 

(1) The ADI technique, although not unconditionally stable as a linear Von Neumann 
analysis indicates, does allow time steps many times larger than most consistent expli- 
cit methods. 

(2) Choosing a criterion for convergence which is too strict may result in no 
steady-state solution being reached. 
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TABLE I.- RESULTS FOR R = 100 


N ' 

Vorticity at midpoint 
of moving wall 

Maximum 
stream function 

15 

8.9160 

-0.08742 

17 

8.4646 

-.09098 

a 17 

7.3756 

-.09867 

33 

6.9919 

-.10038 

57 

6.6960 

-.10128 


a Divergence form. 




TABLE n.- CALCULATED VORTICITY AND STREAM FUNCTION FOR SQUARE CAVITY WITH R * 100 AND EQUALLY SPACED 57 x 57 GRID 


(a) Stream function 











TABLE H.- CALCULATED VORTICITY AND STREAM FUNCTION FOR SQUARE CAVITY WITH R » 100 AND EQUALLY SPACED 57 X 57 GRID - Concluded 


(b) Vorticity 
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-8.093 x 

10 J 2 

1.470 

X 

10-1 

3.443 

x-10-1 

5.181 

X 

10-1 

6.401 x 

10-1 

6.516 x 

10-1 

4.793 

X 

10-! 

7.067 

X 

10-2 

-5.260 

X 

10-1 

-1.091 

X 

10+0 

-1.325 

X 

10+0 

-1.112 

X 

io+o 

.3571 

-1.004 x 

10+0 

-6.322 

X 

10*1 

-3.435 

X 

lO'l 

-1.462 x 

10-1 

-3.450 

X 

10-3 

1.126 

x 10- 1 

2.018 

X 

10-1 

2.388 x 

10-! 

1.865 x 

10-! 

1.356 

X 

lO" 2 

-2.761 

X 

10-1 

-6.029 

X 

10-1 

-8.463 

X 

10-1 

-8.533 

x 10-1 

-6.003 

X 

10-1 

.2857 

-6.233 x 

10- 1 

-4.492 

X 

10- 1 

-3.008 

X 

10' 1 

-1.995 x 

10*1 

-1.298 

X 

10-1 

-7.603 

x 10- 1 

-4.025 

X 

10-2 

-4.022 x 

10-2 

-9.653 x 

10-2 

-2.169 

X 

10- 1 

-3.807 

X 

10- 1 

-5.323 

X 

IQ- 1 

-5.961 

X 

10- 1 

-5.159 

x 10- 1 

-2.831 

X 

10- 1 

.2143 

-3.115 x 

10-1 

-2.927 

X 

10*1 

-2.572 

X 

10-1 

-2.367 x 

10-1 

-2.292 

X 

10-1 

-2.241 

x 10- 1 

-2.203 

X 

10-1 

-2.275 X 

10-1 

-2.572 x 

10- 1 

-3.114 

X 

10-1 

-3.745 

X 

10-1 

-4.146 

X 

io-i 

-3.950 

X 

10-1 

-2.900 

x 10- 1 

-9.135 

X 

10-2 

.14 28 

-9.455 x 

10-2 

-1.711 

X 

10' 1 

-2.124 

X 

lO'l 

-2.567 x 

10-1 

-3.051 

X 

10-1 

-3.453 

x 10-' 

-3.685 

X 

10-1 

-3.755 x 

10-1 

-3.722 x 

10- 1 

-3.632 

X 

10-1 

-3.464 

X 

10-1 

-3.127 

X 

io-i 

-2.504 

X 

10-1 

-1.481 

x 10- 1 

1.013 

X 

lO- 2 

.0714 

-6.710 x 

10-3 

-8.313 

X 

10-2 

-1.608 

X 

10-1 

-2.579 x 

10-1 

-3.635 

X 

10-1 

-4.559 

x lO" 1 

-5.163 

X 

10-1 

-5.328 x 

10- 1 

-5.029 x 

10- 1 

-4.331 

X 

10-1 

-3.375 

X 

10-1 

-2.339 

X 

10-1 

-1.339 

X 

10-1 

-5.867 

x 10-2 

2.978 

X 

lO" 2 

0 

0 



6.610 

X 

10-3 

-8.180 

X 

10-2 

-2.378 x 

10- 1 

-4.141 

X 

10-1 

-5.763 

x lO -1 

-6.946 

X 

10- 1 

-7.403 x 

10- 1 

-6.929 x 

10- 1 

-5.532 

X 

10-1 

-3.528 

X 

lO " 1 

-1.498 

X 

10-1 

-8.960 

X 

lO " 3 

2.866 

x 10-2 

0 






u through point of maximum \p with - 


15 x 15 grid, 
nondivergence form 


57 x 57 grid, 
nondivergence form 


17 x 17 grid, 
divergence form 


.0625 

.0714 

.1250 

.1428 

.1875 

.2143 

.2500 

.2857 

.3125 

.3571 

.3750 

.4286 

.4375 

.5000 

.5625 

.5714 

.6250 

.6429 

.6875 

.7143 

.7500 

.7857 

.8125 

.8571 

.8750 

.9286 

.9375 


-2.612 x IQ' 2 


-5.018 x 10-2 


-7.583 x IQ' 2 


-1.054 x lO' 1 


-1.388 X 10-1 


-1.724 x 10-1 


-1.968 x 10-1 


-1.980 x 10-1 


-1.616 x 10-1 


-7.891 x lO- 2 


5.022 x lO' 2 


2.445 x 10-1 


5.470 x 10-1 


-3.556 x lO' 2 


-6.774 x 10-2 


-1.019 x 10-1 


-1.409 x lO'l 


-1.834 x 10-1 


-2.215 x 10-1 


-2.404 x 10-1 


-2.233 x 10-1 


-1.606 x 10-1 


-5.416 x 10-2 


9.200 x lO' 2 


2.929 x 10-1 


5.905 x 10-1 


-3.509 x lO* 2 


-6.513 x lO' 2 


-9.419 x lO' 2 


-1.245 x 10-1 


-1.563 x 10 "1 


-1.868 X 10-1 


-2.107 x 10-1 
-2.198 X 10-1 
-2.057 x 10* 1 


-1.622 x 10-1 


-8.745 x lO' 2 


1.782 x lO' 2 


1.578 x 10-1 


3.519 x 10-1 


6.315 x 10-1 


1.0000 


1.0 


1.0 


1.0 
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6. COMPARATIVE STUDY OF TWO NUMERICAL TECHNIQUES FOR 
THE SOLUTION OF VISCOUS FLOW IN A DRIVEN CAVITY 

Robert E. Smith, Jr., and Amy Kidd* 

Langley Research Center 

SUMMARY 

Numerical solutions have been obtained for the steady two-dimensional flow of a 
viscous incompressible fluid in a rectangular cavity. The Reynolds number has been 
varied from 100 to 5000. The conservative stream -function — vorticity form of the two- 
dimensional incompressible Navier -Stokes equations was solved. Two techniques 
(alternating -direction implicit and hopscotch) have been applied to the vorticity equation. 
The stream -function equation has been solved by successive overrelaxation and by the 
Buneman direct method. The investigation assesses which combination of these tech- 
niques, applied to the stream -function —vorticity form of the Navier -Stokes equations, is 
the most computationally efficient. Width-depth ratios of 1, 2, and 4 have been used in 
this numerical experiment. Relative efficiency has been based on the central processing 
time of each technique run on the same computer. The results indicate that the 
alternating-direction implicit technique for solving the vorticity equation and the 
Buneman direct method for solving the stream -function equation yield the most efficient 
combination of the techniques studied. In addition, solutions have been obtained at 
Reynolds numbers of 100, 1000, and 5000 with a second-order accurate approximation. 

INTRODUCTION 

Numerical solutions of the stream -function — vorticity form of the Navier -Stokes 
equations describing two-dimensional flow have been the subject of many investigations. 

In particular, the flow properties of a viscous incompressible fluid within a rectangular 
cavity whose upper surface is translating with uniform velocity have been established for 
low Reynolds numbers. (See paper no. 1 by Rubin and Harris.) This study has attempted 
to increase understanding of optimum solution techniques for solving the two-dimensional 
Navier -Stokes equations at low Reynolds numbers. The optimum solution technique is 
taken to be the one which provides a steady-state solution in minimum computer time 
with the least complex program. It is recognized that optimum in this sense is computer 
dependent. However, the results should hold for any third generation computer. Two 
techniques which appear to be near optimum for solving the vorticity equation are the 
alternating -direction implicit (ADI) and hopscotch (HS) procedures. Reference 1 should 
be referred to for the ADI technique and reference 2 is adequate for the HS technique. 

♦Computer Sciences Corporation, Hampton, Va. 
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For solving the stream -function equation the successive overrelaxation (SOR) technique 
and the Buneman direct method (BDM) are candidates for being optimum solution proce- 
dures. Reference 3 should be consulted for the BDM and reference 2 is adequate for the 
SOR technique. 

The two-dimensional Navier -Stokes equations have been described in conservation 
form and approximate boundary conditions established. A computer program with the 
option of using the ADI or HS technique for the solution of the vorticity equation has been 
written. In addition, the HS technique is applied in two ways. The first approach (HS1) 
is direct application of the difference equations at the even and odd grid points. The 
second approach (HS2) is application of a combined difference equation when the iteration 
plus the grid point is even. The second approach should reduce the number of operations 
by approximately 40 percent. The computer program has the option of solving the 
stream -function equation by either SOR or BDM. The BDM subroutine used in the com- 
puter program is one written by Buzbee. (See ref. 3.) The Buneman algorithm is much 
more complex than the SOR algorithm; however, once the algorithm is written into a 
computer program, it is easily applied. Note that if the domain of interest is not rec- 
tangular, the Buzbee routine and the BDM are not directly applicable. The computer 
program has been written with the ability to vary time step, grid size, Reynolds number, 
and initial conditions. With this program as an experimental tool, many runs varying 
time step, grid size, Reynolds number, and width-depth ratio with each of the combina- 
tion techniques - ADI-SOR, ADI -BDM, HSl-SpR, HS1-BDM, HS2-SOR, and HS2-BDM - 
have been made on a CDC 6400 computer and appropriate statistics recorded. 

SYMBOLS 

N,M dimensions of grid (see fig. 1) 

R Reynolds number 

t time 

At time step 

U uniform velocity at upper surface (see fig. 1) 

u axial velocity 

v- - 

v normal velocity 
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x,y axial and normal coordinate, respectively 

Ax, Ay spatial increments in x- and y-directions 
/3 = Ax/Ay 

£ vorticity 

i// stream function 

a; convergence parameter in SOR approximation 

Superscripts: 

k index of iteration 

n index of time step 

Subscripts: 

i,j index of grid point in x- and y-direction, respectively 

in grid point inside surface 

out grid point outside surface 

s grid point on surface 

t differentiation with respect to t 

x,y differentiation with respect to x or y 

Abbreviations: 

ADI alternating direction implicit 

BDM Buneman direct method 

HS1 hopscotch method 1 
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HS2 


hopscotch method 2 


SOR 


successive overrelaxation 


GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 


The normalized stream -function —vorticity conservation form of the two- 
dimensional incompressible Navier -Stokes equations is: 


?t = -(^y)Sx + (^x)^y + r(?xx + ?yy) 


( 1 ) 


xb + xb = _£ 
XX T *yy s 


( 2 ) 


The boundary conditions on the stream -function equation (eq. (2)) for flow in a 
cavity with the upper surface translating to the left with uniform velocity U = 1 and 
with no flow at the other boundaries are (see fig. 1): 

at upper surface (moving to left), 


\b = _1 


*x = ° 


\p= 0 


at bottom, left, and right surfaces, 


iZ/ =0 


*x = ° 


Ils = 0 


The boundary conditions on the vorticity equation are obtained by central differ- 
encing equation (2), by applying the boundary conditions for the stream function, and by 
enforcing reflection at the boundary. They are: 

at upper surface, 


^out ~ + ^in _ 


(Ay)' 


= -S* 


*s = ° 


^out ~ ^in _ 
2 Ay 
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or combining these conditions results in 


2 ^ in - 2 Ay 


(Ay) 


2 


at bottom surface, 


^out - IV^in _ 


(Ay)' 


= 


^ = 0 


^in ^out 


or 



s 


at left and right surfaces, 


^out 


- 2 ^s + ^in _ • - 
(Ax) 2 


s 


^s = ° 

^in = ^out 


or 


2^; 


in 


(Ax) 2 



FINITE-DIFFERENCE APPROXIMATIONS 


Because of the complexity of the BDM, an explanation of this method is not given. 
Details of the method are given by Buzbee, Golub, and Nielson in reference 3. 
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ADI Approximation for Vorticity Equation 
The ADI algorithm is applied to the vorticity equation (eq. (1)) in the following two 

steps: 

Step 1 


n+1/2 Vx * 




f" 4 , 1 / 2 + (1 + 2X x )C n+1 / 2 + 

^ X ^1,J 




j»n+l/2 

M+1J 


■My-* + ^ 




y > n 

C U+1 


(3) 


Step 2 


Wyli^y - V 


«y -i + t 1 + 2 x vK] L + 


n+l 




C 


n+l 


i+l,j 


" +1/2 y„ + V 


TO 


X T "X 


cf-if + (i - 2> xX"i 1/2 


+ 



1/2. 

,j 





(4) 


where 


At 

4 Ax 


_ At 
4 Ay 


_ At 
2R(Ax) 2 


A - At 
y 2R(Ay) 2 


and 


Ax = X[ - Xj_i 


Ay = yj-yj_i At = t n - t n_1 


Steps 1 and 2 form tridiagonal systems of linear algebraic equations. Solving 
equations (3) and (4) in succession yields one time step. 
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Hopscotch Approximation for Vorticity Equation 
The following three equations are used in the hopscotch techniques: 


j»n+l 


N 




y> n 




3 1 J 


i+M + 




V 


'",M + (>- 2x i - 


( 5 ) 


'y 1 


'1-1J 


J»n+1 , 

c i-i,i + 




j.n+1 
■>+ 1 j 




<p n 

? U-1 


Wi,j!i v y + x y 


rP . , + r!*.l 
s i,J+l ^1,3/ 


( 6 ) 


and 


'u 1 = 2 ?u - ?fj‘ 


(7) 


where 


\ ~ At 

A. x 


R(Ax)‘ 


, _ At 


x y = 


R(Ay)‘ 




- At 
2 Ax 


y- = ^ 

f V 


y 2 Ay 


The HS technique is used in two ways to approximate the vorticity equation: 

(1) In the first method (HS1), the difference equations are directly applied in two 
sweeps. First, equation (5) is computed at grid points for which i + j + n is even. 
Then, equation (6) is computed at points for which i + j + n is odd. 

(2) The second approach (HS2) is a three -level method which results from a com- 
bination of equations (5) and (6). For n = 2, the equations are applied as for HS1. For 
n = 3, equation (7) is computed at points for which i + j + n is even and equation (6) is 
applied at points for which i + j + n is odd. 
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The hopscotch technique is conditionally stable for the vorticity equation according 
to the Courant -Friedrichs -Lewy (CFL) condition; that is, U(At/Ax) ^ 1. For HS1, 

21 arithmetic operations are necessary for computing C at time level n + 1. Dimen- 
sional computer variables are stored and retrieved from memory at two time levels, n 
and n + 1. With HS2, 43 percent less arithmetic operations are necessary. However, 
HS2 requires 33 percent more storage than HS1. Since HS2 requires storage of two- 
dimensional arrays at three time levels, 43 percent less arithmetic operations does not 
necessarily imply 43 percent increase in computational efficiency. 

SOR Approximation for Stream -Function Equation 

The following equation is used in the SOR routine to approximate the stream 
function: 

*k + l . *k . + + - (ta)2 ?y - 2(l ♦ 02)^' 

2\1 + /3 ) L 

where the convergence parameter o> is 1.67351. 

ANALYSIS AND RESULTS 

An analysis of the techniques discussed in the preceding section based on an exper- 
imental computer program designed to execute the ADI-SOR, ADI-BDM, HS1-SOR, 
HS1-BDM, HS2-SOR, or HS2-BDM combination has been performed. The results of com- 
puter runs made with this program on a CDC 6400 computer are summarized in table I. 
The Reynolds number has been varied from 100 to 5000 in cavities with width-depth of 1, 
2, and 4. The basic grid size is 1/16. More accurate solutions were obtained with grid 
sizes of 1/32 and 1/64 for a Reynolds number of 100 and a width-depth ratio of 1. In 
order to save computer time with the 1/64 grid, initial conditions with the converged 
1/32 grid were used; Linear interpolation was used to obtain initial conditions at the 
midpoints. The 1/64 grid required considerably more computer resources than any of 
the others. The initial conditions for all runs other than those with the 1/64 grid are 
i// = 0.1 and £ = 0. 

For the solution of flow in a cavity, the ADI-BDM combination was found to be the 
most computationally efficient of the techniques studied. The memory requirements are 
comparable with the HS1-BDM combination which is approximately 20 percent slower. 

In both the ADI-BDM and HS1-BDM combinations, only two levels of information must be 
stored simultaneously. The advantage of the ADI-BDM combination is due to the larger 
time steps that can be taken. The HS1-BDM and HS2-BDM combinations require fewer 
arithemetic operations per time step; however, hopscotch is stability limited by the 
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CFL condition (U(At/Ax) = l). A linearized stability analysis indicates the ADI method 
is unconditionally stable. In practice, arbitrarily large time steps could not be used. It 
was necessary to make runs with a number of time steps to determine a satisfactorily 
large step which produced a stable solution for a particular Reynolds number and grid 
size. 

The BDM is probably the fastest technique presently known for solving a Poisson 
equation over a rectangular domain. The particular routine of Buzbee (ref. 3) is partially 
written in COMPASS (machine language) to add additional speed to the routine. Since 
this routine yielded a 1/3 decrease in central processing time, it was used exclusively 
for more time-consuming runs. 

The HS1 method was relatively easy to program and its storage requirements sat- 
isfactory (two simultaneous levels of information). It was not necessary to run this 
technique with At = Ax since the velocity falls off from 1 to 0.6 (R = 100) in one grid 
step (Ax = Ay = 1/16). However, when At > Ax, some oscillation would occur near the 
moving boundary. 

The HS2-BDM combination was the second most computationally efficient on the 
basis of these computer experiments. It was about 20 percent slower than ADI -BDM 
combination. According to the number of arithmetic operations, it should have required 
about 40 percent less computational time than HS1-BDM; however, because of the 
requirement of three simultaneous time levels of information, there are added machine 
operations associated with storing and retrieving this information. Table I shows the 
results of all the computer runs in order of ascending central processing time. Because 
of its explicit nature, it is believed that the hopscotch algorithm has a high potential for 
the STAR computer. It is 20 percent slower than the ADI method, but the machine effi- 
ciency of an explicit (vectorizable) technique has the potential of rendering HS1 superior 
to ADI when an inconsistent transient solution is acceptable. 

REMARKS ON ACCURACY 

Solutions calculated with the ADI technique of the stream -function — vorticity form 
of the Navier -Stokes equations are given in figures 2, 3, and 4 for Reynolds numbers of 
100, 1000, and 5000, respectively. The results for a Reynolds number of 100 (fig. 2) 
agree well with Burggraf (ref. 4) and Bozeman (ref. 5). For a Reynolds number of 1000, 
the results for the stream function obtained in this study (figs. 3(a) and 3(b)) have the 
character of those obtained by Bozeman. Bozeman used a grid size of 1/50; however, 
his differencing of the convective terms of the vorticity equation for a Reynolds number 
of 1000 is first-order accurate, whereas the present solution is based on a central - 
difference approximation of the convective terms and is therefore second-order accurate. 
The contour plots of the vorticity solutions for Reynolds numbers of 1000 and 5000 
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(figs. 3(c) and 4(c)) are meaningless near the center of the cavity since the surface is 
nearly flat there. Probably the perspective plots (figs. 3(d) and 4(d)) are more mean- 
ingful. They show how. the vorticity changes rapidly at the boundaries and is flat in the 
middle for high Reynolds numbers. The 1/16 grid size is too large for good accuracy 
at Reynolds numbers of 1000 and 5000 even though the solutions obtained show the 
expected character. A point of interest is that with the ADI procedure, which used cen- 
tral differencing of the convective terms in the vorticity equation, converged solutions 
could not be obtained at higher Reynolds numbers than Bozeman used. Bozeman was 
unable to obtain solutions with either central or upwind differencing for R > 1000. 

CONCLUDING REMARKS 

The primary conclusion of this study is that a combination of the alternating - 
direction implicit technique and the Buneman direct method is a rapid procedure for 
solving the Navier -Stokes stream -function and vorticity equations in a rectangular 
region. The combination of the hopscotch and Buneman direct methods is only 20 per- 
cent slower at low Reynolds numbers. The hopscotch technique. was applied in two ways: 
(1) direct application of the difference equations at even and odd grid points and (2) appli- 
cation of a combined difference equation at points for which the iteration plus the grid 
point was even. With both hopscotch approaches combined with the Buneman direct 
method, the first is preferable to the second, because it is easier to program, requires 
less storage, and is only slightly slower. The first hopscotch approach appears to have 
potential for flow computations on vector machines such as the STAR computer. 

It was possible to obtain converged solutions at Reynolds numbers of 1000 and 
above with the alternating -direction implicit technique which used central -difference 
approximations for the convective terms in the vorticity equation. The solutions obtained 
at higher Reynolds numbers are not necessarily accurate, but they do show the character 
of high vorticity gradients near the moving boundary and low vorticity gradients near the 
center. 
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TABLE I.- RESULTS OF ANALYSIS OF TECHNIQUES ON CDC 6400 COMPUTER 


[Convergence criterion was 0.0001] 


Reynolds 

number 

Width-depth ratio 
of cavity 

Techniques used - 

Grid size, 
Ax = Ay 

Time 

step 

Central 
processing 
time, sec 

Number 

of 

iterations 

Computer 

time 

ratio 

(a) 

For vorticity 

For stream function 

100 

1 

ADI 

BDM 

1/16 

0.2 

56.85 

119 

1 

100 

1 

HS2 

BDM 

1/16 

.063 

63.8 

337 

1.12 

100 

1 

HS1 

BDM 

1/16 

.063 

65.4 

339 

1.14 

100 

1 

ADI 

SOR 

1/16 

.2 

90.17 

110 

1.586 

100 

1 

HS2 

SOR 

1/16 

.063 

117.8 

325 

2.07 

100 

1 

HS1 

SOR 

1/16 

.063 

122.9 

328 

2.16 

100 

2 

ADI 

BDM 

1/16 

.2 

260.9 

136 

4.58 

100 

4 

ADI 

BDM 

1/16 

.2 

588.8 

139 

10.35 

100 

1 

ADI 

BDM 

1/32 

.075 

478.7 

322 

8.42 

100 

1 

ADI 

BDM 

1/64 

.0097 

16 783.0 

2561 

295.21 

1000 

1 

ADI 

BDM 

1/16 

.2 

198.03 

693 

3.48 

5000 

1 

ADI 

BDM 

1/16 

.2 

2 106.5 

3469 

37.05 


a Ratio of central processing time to central processing time of first entry. 
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(e) Velocity field. 
Figure 3.- Concluded. 





X = 1 

x = 0 


(a) Contour plot of stream function. 
Figure 4.- ADI-BDM solutions for R = 5000. 
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7. CALCULATION OF STEADY VISCOUS FLOW IN A SQUARE DRIVEN 


CAVITY BY THE ARTIFICIAL COMPRESSIBILITY METHOD 

John T. Suttles 
Langley Research Center 

SUMMARY 

Numerical solutions of steady viscous flow in a two-dimensional cavity have been 
obtained by the method of artificial compressibility. The method, originally presented 
by A. J. Chorin, was programed for a computer solution and checked against a simple 
channel flow problem which has an analytical solution. The comparison indicated that 
the method was correctly programed. Results for a square cavity at a Reynolds num- 
ber of 100 were then obtained and compared with published results. This comparison 
showed that the method produces good results for velocities but that the computed pres- 
sures are unreliable. 


INTRODUCTION 

The literature, in recent years, includes numerous studies of different numerical 
methods for obtaining finite-difference computer solutions to fluid mechanics problems. 
(See ref. 1, for examples.) Consequently, a study was undertaken to compare the ability 
of several reported methods to compute the motions of a viscous incompressible fluid in 
a two-dimensional square cavity which has one of its walls undergoing uniform transla- 
tion. (See paper no. 1 by Rubin and Harris.) This problem was selected because it rep- 
resents a simple example of a steady flow involving closed streamlines and it has been 
studied extensively in the literature. (See refs. 2 to 6.) The present paper is one of a 
series of papers dealing with this problem and reports the results of application of the 
method of artificial compressibility (ref. 7). 

The method of artificial compressibility was proposed by A. J. Chorin (ref. 7) as 
an efficient approach to solving the steady-state problem of viscous incompressible flow. 
The principle of this method is to solve the steady incompressible problem by using an 
asymptotic time solution of the equations of motion. These are modified to contain an 
artificial compressibility that is designed to vanish as steady state is approached. 
Through a stability analysis and numerical experiments, Chorin showed that the approach 
involves two parameters (time step At and artificial compressibility 6) which may be 
adjusted within the constraints of the stability criteria to obtain a converged solution rap- 
idly. In reference 7, a simple channel flow was used as a test problem, and the method 
was applied to the solution of thermal convection in a fluid layer heated from below. 
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Because of the analytical complexities involved, Chorin relied heavily on his numerical 
results to show the stability, convergence (to a steady state in time and to a unique solu- 
tion as the grid size is reduced), and efficiency of the method. 

In the present paper, the problem of steady viscous flow in a square cavity is form- 
ulated in primitive variables (as opposed to a stream -function — vorticity formulation), 
and the method of artificial compressibility is applied to obtain the solution. Therefore, 
these results provide further numerical evidence for evaluating the method and a con- 
sistent basis for comparing the accuracy and efficiency of this method with other 
approaches. 

SYMBOLS 


Cp pressure coefficient 

c artificial speed of sound 

l characteristic length 

A l = Ax = Ay 

M artificial Mach number 

N number of space dimensions 

p pressure 

/ n o\l/2 

q magnitude of velocity vector, (u* + v 6 ) 

R Reynolds number 

t time 

At time step 

U characteristic velocity, velocity of translating wall of cavity 

u axial velocity 
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V 


normal velocity 


x,y axial and normal coordinate, respectively 

Ax, Ay spatial increments in x- and y-directions 

6 artificial compressibility 

u kinematic viscosity 

p density 

Superscript: 

n index of time step 

Subscripts: 
c center 

i,j index of grid point in x- and y-direction, respectively 

max maximum value 

ref reference 

A bar over a symbol denotes a dimensional quantity. 

METHOD 

Governing Equations 

The equations of motion of an incompressible viscous fluid in two-dimensional 
rectangular coordinates are: 

Continuity 

91 + 91=0 ( 1 ) 

ax ay 
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x -momentum 


9M + u-^ + v^ = -1^1+ v($-± + 
at 3x 9y p ax ^2 9 y2 


( 2 ) 


y-momentum 

a| + n-5l + v^ = -1 ^ 

at ax ay p ay [ dx 2 dy 2 j 


(3) 


where bars over symbols denote dimensional quantities. These equations are nondimen - 
sionalized by defining: 


u 






P 


= = P 
pv U 



where l is a dimensional characteristic length and U is a dimensional characteristic 
velocity. Applying these relations in equations (1) to (3) and using the continuity equa- 
tion to write the momentum equation in divergence form result in 


+ o 

ax ay 


au . p/au 2 , auv\ _ 9 P , /a 2 u , a 2 u\ 
at lax ay/ "'ax (^2 dy 2) 


av , p/auv , av 2 l _ ap _ /a2 v _ g2 v ^ 
at lax ay/ “ay [ dx 2 + dy 2j 


where the Reynolds number R is defined by 



v 


(4) 

(5) 

( 6 ) 


(7) 
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To apply the method of artificial compressibility, the continuity equation (eq. (4)) 
is rewritten as 


9 P 
3t 




( 8 ) 


where p is the artificial density and t is an auxiliary time which is analogous to real 
time in a compressible flow problem. An artificial equation of state, 



is introduced with 6 being the artificial compressibility. Equation (8) may be solved 
with equations (5) and (6) to produce a pseudotransient solution. Note that in such a 
solution the time in equations (5) and (6) becomes the auxiliary time t . The solution 
thus obtained is meaningful only if a steady state is attained and the time derivatives in 
equations (5), (6), and (8) vanish so that p and t no longer influence the results. 
Thus, this approach is appropriate when a steady-state solution is desired. 


Finite -Difference Equations 

In applying the method of artificial compressibility, the governing equations may be 
written according to any of the known difference schemes. The schemes selected by 
Chorin (ref. 7) were used in the work reported herein. These are the leapfrog difference 
scheme, wherein central differencing is applied to first-order time and space derivatives, 
and the DuFort-Frankel method for the second-order space derivatives (viscous terms). 
The DuFort-Frankel approach treats the derivatives as follows: 


c% nH . till 

a 2 u _ i+l,J i-l,j 


- u 


n+1 _ u n-l 




i,l 


ax* 


(Ax) 2 


where 

u(x,y,t) s uf ; j = u(i Ax,j Ay,n At) 


and Ax, Ay, and At are the space and time increments. This scheme is discussed in 
reference 1, in which it is shown that use of the DuFort-Frankel scheme results in a dis- 
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torted time scaling. However, since the method of artificial compressibility seeks only 
the steady-state solution, this difficulty with the DuFort-Frankel approach is of no 
consequence. 

With the above schemes, the finite -difference form of equations (8), (5), and (6) 
becomes 


m+1 

i,j 


- pP" 1 

"i ,] 


2 Af 


„n 

u i+l,j 


- u 


i-l,j 


2 Ax 


,n 


i>j+l 


2 Ay 


n 

i,3-l 


u 


n+1 „n-l 


i,l 


- u 


i,3 _ 


= -R 


2 At 


Km) 2 - ( u "-i,i) 2 

2 Ax 


R 


u i!j+l vI1 


U+l 


- U P- lV P. 

i,J -1 i,J- 


2 Ay 


} i + l,j -Pi-l,J 

2 Ax 


uP , • + uP , . - uP+ x - uPtI uP . , + uP . , - uPt 1 - uPt 1 

. t+1,3 1-1,1 i,3 i,3 . M+l i,l-l i,J i,J 

« A ' 


(Ax)' 


(Ay)' 


V P+1 - yP' 1 
1,1 1,3 = 

2 At 


- u i I -i,j v "-i,i Ki+i) 2 - Ki-i) 2 i ff,j+i 

2 Ax * 2 Ay ' 6 2 Ay 


,n : : x „n 


r n+l TJ n-l ,,n 


,n+l ,,n-l 


v. t . + v- 1 • - V* *. - V- • V. • 1 4- v. . .. - v. . _ V. . 

i+l,] i-l,J i,3 i,3 1,3+1 1,3-1 i,J 1,3 


(Ax)' 


(Ay)' 


where uP^ = u(i Ax,j Ay,n At ), vf^ = v(i Ax,j Ay,n Af), and pP^ = p(i Ax,j Ay,n Af) and 
equation (9) has been used to eliminate pressure from the system of equations. Assum- 
ing equal space increments Ax = Ay = A l and solving for the artificial density and 
velocities explicitly result in 



(10) 


£(pP 
6V h 


i+l,j 



, 2^/ u n , 

aA i+1,3 


+ u y+i + 


Hi) 
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( 12 ) 


v-H* = avP;: 1 - R/3 
l jJ 




+ M/ V P , . + v? , • 

M\ i+l, J i-l,3 


+ v y+i + 


vP. ,\ 


where 


1 -4-^L 

(AZ ) 2 

1 + 4 “^2 

(M ) 2 


(13) 


and 


0 = 


Af 

A l 


1+4 


At 

(AZ ) 2 


(14) 


The boundary conditions consist of prescribed velocities on all boundaries; however, 
the artificial density given by equation (8) must be updated on the boundaries. In refer- 
ence 7, Chorin used forward or backward space differences at the boundaries to update 
the artificial density in the continuity equation (eq. (8)). For example, at a point on the 
y = 0 (j = 1) boundary 


^n +1 _ ^n -1 At /.,n n n \ o At / v r i v n A 

p i,l ~ p i,l ' azTI+M " " 2 a 7\ b 2 ' Vi ^) 


(15) 


Points on other boundaries were treated similarly. This method of treating the boundary 
points was initially adopted in this work also. The approach was found to work quite 
satisfactorily for the duplication of the simple channel flow problem presented by Chorin; 
however, it was found to produce poor results near the boundaries for the cavity problem. 
An attempt to rectify this problem is described in the section entitled "Results and 
Discussion." 


Stability 

The stability conditions for the system of finite -difference equations (eqs. (10) to 
(12)) have been discussed by Chorin (ref. 7). The main points of that discussion are given 
in this section. 
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The artificial equation of state, equation (9), indicates the presence of an artificial 
sound speed 


c = 


R6 1 / 2 


(16) 


where the Reynolds number has been introduced to allow c to be nondimensionalized in 
the same manner as the other velocities. To ensure the stability of the system of equa- 
tions, the artificial Mach number based on this sound speed must be less than unity. 

The required condition is 


“max - 5sp = R« 1/2 q ma x < 1 


(17) 


where q is the magnitude of the velocity vector; that is 



+ v 


1/2 


If this Mach-number condition is satisfied and if the boundary conditions consist of pre- 
scribed velocities, Chorin showed that the remaining stability requirement is satisfied 

by 


At S 



(18) 


where N is the number of space dimensions and ^min Axj^ is the minimum spatial 

increment. Thus the method contains two parameters, 6 and At, which may be 
adjusted within the stability criteria given by equations (17) and (18), respectively, to 
optimize the time required for obtaining a converged solution. Chorin states that the 
optimal value of 6 must be determined by preliminary test computations and that his 
experience indicated that this value is not sharply defined. Therefore, no attempt to 
determine optimal values was made in this work. It was found, however, that a maxi- 
mum Mach number of 0.5 produced satisfactory results in most instances. 

RESULTS AND DISCUSSION 

Simple Channel Problem 

In reference 7, Chorin presented results of calculations for a simple two- 
dimensional channel flow problem which has a closed form solution, that is, a parabolic 
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axial velocity profile across the channel and a linear decay of pressure going down the 
channel. Calculations for this test problem have also been made in this work to gain con- 
fidence that the method has been correctly programed. Both a horizontal (channel axis 
parallel to x-direction) and a vertical (channel axis parallel to y-direction) channel flow, 
shown in figure 1, were computed. In these calculations, the boundary velocities were 
specified to be zero except for the axial velocities across the channel ends which were 
specified to have the correct parabolic distribution. The interior velocities and densities 
were initially set at zero, and the velocities and densities in the interior of the channel 
were computed. The calculations were done for a unit Reynolds number (R = 1) with 
19 grid points in each coordinate direction and for 6 = 0.000 32 and At = 0.000 596 28 
as was done in reference 7. The computations for both channels converged to the cor- 
rect parabolic axial velocity distributions (i.e., all velocities within 10“^ of known solu- 
tion) and to the correct linear axial pressure gradient throughout the channel in approxi- 
mately 500 iterations. This rate of convergence is in good agreement with Chorin's 
reported experience. It was, therefore, concluded that the method has been properly 
programed. 


Square Cavity Problem 

The square cavity problem considered in this work is illustrated in figure 2. The 
steady viscous flow in the cavity has been computed by Chorin's artificial compressibil- 
ity method. Preliminary computations for this problem indicated that the velocity 
results near the boundaries were poor. The difficulty was found to be due to the use of 
conditions like equation (15) to update the artificial densities along the boundaries. This 
problem was eliminated by applying the governing equations at the boundary to derive a 
relationship for the pressure gradient normal to the boundary (and from the artificial 
equation of state, the artificial density gradient also) in terms of the normal velocity 
gradient. For example, along the boundary where y = u = v = 0, and hence 
3u/9t = 9v/9t = 9u/3x = 3v/9x = 0, the relation is found from equation (6) to be 



(19) 


Taylor's series expansions about the boundary point for the pressure and velocity gradi- 
ent one grid point from the boundary result in 


p(i,2) = p(ijl) + Ay 

\ a y/i,i 
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Equation (19) can then be used along with these relations to obtain 
P(M)=P(i,2)-(H). j2 + (H) M 

Finally, with a central -difference expression for the velocity gradient, the boundary pres 
sure is given by ✓ 

p(i,l)=p(i,2) (20) 


when it is noted that 
v(i,l) = 0 

from the boundary condition and that 

= 0 
i,l 

from the boundary condition applied to the continuity equation. Thus, the appropriate 
equation for updating the artificial density on the y = 0 (j = 1) boundary is found from 
equation (20) and the artificial equation of state (eq. (9)) to be 

p(i,l)=p(i,2) (21) 



Analogous equations for the remaining boundaries are: 
for the x = 0 (i = 1) boundary, 

3(1, i) - 3(2, j) - 

2 A l 


( 22 ) 
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for the y = 1 (j = j max ) boundary, 


P(Umax) “ P^>imax - 1) + 


5 v(i,j 


max 
2 A l 


2 ) 


(23) 


for the x = 1 (i = i max ) boundary, 

P(imax>j) = P(^max - l>j) + — (24) 

2 A l 

The use of equations (21) to (24) made a considerable improvement in the velocity results 
near the boundaries. 

Calculations for a square cavity were made for a Reynolds number of 100 with 
15 x 15, 29 x 29, and 43 x 43 grid points. Table I lists the stability parameters, 6 and 
At, the number of iterations, required storage, and run time for the calculations. 

Although a large number of iterations were required (1500 to 3000), the computations 
only took computational times of approximately 1/2, 3, and 7 minutes for 15 x 15, 29 x 29, 
and 43 x 43 grid points, respectively. The results of these calculations are presented 
and compared with other results in figures 3 to 7. 

In figure 3, the velocity vectors are shown for the steady-state solution by the 
method of artificial compressibility with 15 x 15 grid points. The flow direction vectors, 
which are obtained by normalizing each velocity vector to a unit magnitude, are shown in 
figure 4. These results illustrate the features of the vortical flow which occurs for the 
cavity problem. The features shown are in good agreement with the streamline patterns 
shown in reference 3 for this problem. The anomalous flow direction vector shown in 
the lower right corner of figure 4 is probably associated with a weak corner vortex flow 
which is described in references 3 and 5. 

The velocity and pressure profiles through the vortex center (i.e., the nearest node 
point in the grid to the vortex center) for the present solutions are compared with the 
results of reference 3 in figures 5 and 6. The velocity results in figure 5 show that the 
present solutions for the various grid spacings vary only slightly and that the results for 
15x15 grid points appear to be adequate. It is also seen that the present calculations of 
velocity are in good agreement with those of reference 3 which were obtained by a 
stream -function — vorticity method. The vortex center locations are apparent in these 
graphs and are compared with those of reference 3 in table II where it seems that the 
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values are in good agreement, that is, within about 4 percent. On the other hand, the 
pressure profiles through the vortex center shown in figure 6 are in substantial disagree- 
ment. Of particular concern is the fact that the pressure profiles for the various grid 
spacings of the artificial compressibility method are not smooth and are in disagreement 
with each other. This makes the pressure results of this method highly questionable. 

In figure 7, the pressures along the boundaries are compared qualitatively with the 
results in reference 6. Note that the Reynolds numbers are not the same, that is, 100 
compared with 64. Although there are relatively high frequency oscillations in the pres- 
ent results, the qualitative comparisons show good agreement with reference 6. 

Thus, these results indicate that the method of artificial compressibility applied to 
the square cavity problem produces good results for velocity, but the pressure results 
are unreliable. 


CONCLUDING REMARKS 

The method of artificial compressibility has been applied to obtain numerical solu- 
tions of the steady viscous flow in a two-dimensional square cavity. The method, pro- 
posed by A. J. Chorin, was programed and checked against a simple channel flow problem 
which has an analytic solution. The comparison indicates that the method was correctly 
programed. Results for a square cavity with a Reynolds number of 100 were obtained 
and compared with published solutions for this problem. These comparisons indicated 
that good velocity results are obtained with the method, but that the pressure results are 
unreliable. 
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TABLE I.- STABILITY AND COMPUTATIONAL PARAMETERS 
FOR SQUARE CAVITY CALCULATIONS 


Number of 
grid points 

6 

At 

Number of 
iterations 

Execution 
time, a sec 

Execution 

storage a 

15 x 15 

0.000 025 

0.000 21 

1500 

29 

20624g 

29 x 29 

.000 003 2 

.000 038 33 

3000 

205 

24314g 

43 x 43 

.000 025 

.000 0714 



3000 

419 

32234g 


a For CDC 6600 computer system which has a minimum storage of 40000g for 
compilation. 


TABLE II.- COMPARISON OF VORTEX CENTER LOCATIONS 


Vortex 

Results 

center at - 

Present 

Mills (ref. 3) 

X 

0.61 


0.63 

y 

.73 


.75 
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0*xSl, y=0 


u - 0, v = 0 
(a) Horizontal channel. 



0 - 0 £ x S 1, y = 0 


u = 0, v = Vx( 1 - x ) 

(b) Vertical channel. 

Figure 1. - Geometry and boundary conditions for simple channel flow problems. 
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Figure 3.- Velocity vectors from solution of square cavity problem 
with R = 100 and 15 x 15 grid points. 


1 t 
\ 1 


1 t 
1 *S 
1 \ 
1 \ 
t \ 
\ \ 
\ \ 


It t T/s 
t I H\ \ \ i i 


\ \ \ ^ 
\ 

\ \\K^ 


S / J l 

s / / / 
^ / / / 
/ / 

^ / / 

^ / J I 
^ s / / 


Figure 4.- Flow direction vectors from solution of square cavity problem 
with R = 100 and 15 x 15 grid points. 
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Figure 6.- Profiles of pressure coefficient through vortex center from present calculations and from reference 3. 
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8. FINITE -DIFFERENCE SOLUTION FOR THE INCOMPRESSIBLE 
DRIVEN CAVITY FLOW PROBLEM 

Ernest V. Zoby 
Langley Research Center 

SUMMARY 

An alternating -direction implicit (ADI) finite -difference procedure has been 
employed to compute the two-dimensional viscous incompressible flow generated within 
a rectangular cavity by a moving surface at the cavity open end. The pressure through 
the cavity was computed by the Spalding SIMPLE procedure. 

Although the present method does not predict axial and normal velocities as high 
as indicated by more exact methods, the values of the velocities through the vortex cen- 
ter are reasonably predicted and the overall velocity vector diagrams appear satisfac- 
tory. The pressure levels are generally in poor agreement with other predictions. 
However, the pressure gradients must be reasonably well predicted because of the accu- 
rate results obtained for the velocities. 

INTRODUCTION 

The solution of the Navier -Stokes equations easily provides a more than adequate 
test for the viability of various numerical techniques. For a rational assessment of 
these methods, a singular flow model and compatible test conditions are naturally desir- 
able. Numerical as well as experimental results (e.g., ref. 1) exist for the two- 
dimensional incompressible flow field generated within a rectangular cavity by a moving 
surface at the open end. (See paper no. 1 by Rubin and Harris for details of driven cav- 
ity.) Consequently, the driven cavity problem was selected as the test model for an 
investigation to compare the accuracy and applicability of several numerical techniques. 

As part of the aforementioned investigation, this paper presents numerical solu- 
tions based on an alternating-direction implicit (ADI) technique for the two-dimensional 
incompressible Navier -Stokes equations. The equations are programed in the rectangu- 
lar Cartesian coordinate system. For an evaluation of the pressure, the Spalding 
SIMPLE procedure (ref. 2) is used. The Spalding method relates a perturbation in pres- 
sure at each time step to the continuity relation so that at steady state the pressure at 
each node point is determined and the continuity equation is satisfied. A discussion of 
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the present method and the results of calculated test cases are presented. Comparison 
with existing results is also included. 

SYMBOLS 


C P 

l 

P 

P 

R 

t 

At 

U 

u 

V 

x,y 

Ax, Ay 

v 

P 


pressure coefficient 

characteristic length (fig. 1) 
pressure 

pressure correction 
Reynolds number 
time 

time step 

velocity across open end of cavity 

axial component of velocity 

normal component of velocity 

axial and normal distance, respectively 

spatial increments in x- and y-directions 

kinematic viscosity 

density 


Subscripts: 

i, j index of node points in x- and y-direction, respectively 

max maximum value 
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Superscripts: 
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( 5 ) 


lX + u 9y + v 9V-. _^P , 1 fd^y a 2 v\ 
3t 3x 9y - 8y R^ ax 2 dy 2J 


and 


3u + 9v _ 0 
ax ay 


( 6 ) 


where R = \Rfv. The boundary conditions on u and v are 

u(0,y) = u(l,y) = u(x,0) = 0 u(x,l) = 1 

v(0,y) = v(l,y) = v(x,0) = v(x,l) = 0 

However, with equations (4) to (6), there is no simple way to compute pressure. 
For this paper, a semi -implicit method of Spalding (ref. 2) has been employed. The 
equation for this method is 


_ (& P £^P\ = M + 9V 
at \ax 2 ay 2 1 ax a y 


( 7 ) 


where the ap/at is a pseudotime term relating to the iterations and p is a pressure 
correction which is combined with the pressure p from the previous time step for a 
new pressure value at the new time step. The correction p continually updates p in 
the time marching procedure until convergence is obtained at steady state. At steady 
state, p is very small, so that the left side of equation (7) equals zero and continuity is 
satisfied. The boundary condition for equation (7) is that the normal pressure gradient 
is zero. 

Equations (4), (5), and (7) are then applied to the square cavity (fig. 1) to determine 
the flow field generated within the cavity. The basic conditions for these calculations 
are that R = 100 and Ax = Ay = 1/14. For initial conditions, u, v, and p (and p) 
were taken to be zero. 

Equations (4), (5), and (7) are solved by an alternating-direction implicit (ADI) 
finite -difference procedure. With the first one -half time step implicit in the y-direction, 
equations (4) and (5) are given as 
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vO • At 

V i,3 , At 

n+1/2 . 

[i+ At 

u n+l/2 . 

vP • At 
1,3 1 At 

4 Ay 2R(Ay) 2 

U • ■ '■* T 

1,3-1 

R(Ay) 2 

i,3 

4 Ay 2R(Ay) 2 


u n+i/2 

i,j+l 


At 

4 Ax 


(P?+I,j ' p ?-l,j) 


ufi At 




*>3 ( At 

^-1,3 + 

i 

At 

4 Ax 2R(Ax) 2 

J. 

R(Ax) 2 


u 


1,3 


At 




2R(Ax) 


2 4 Ax 


„n 

U i+1,3 


( 8 ) 


and 


V? • At 
l >3 + At 


4 Ay 


2R(Ay)‘ 


_ v n+l/2 + 

[l y At 1 

n+1/2 1 

2 V i,j -1 + 

R(Ay) 2 

V i,3 + 


i,3 


At 


4 Ay 2R(Ay) 2 


v *>+ l { 2 

i,j+l 


"4 A Ay( P M+l “ P ?,j -l) + 


i At ' 

*>3 At 

4 Ax 


2R(Ax) 2 


v" .+ 

1-1,3 


1 - 


At 


R(Ax)' 


i,3 


At 


u u At 


2R(Ax) 2 4 Ax 


v p 1 . 
1+1,3 


0) 


For the second one -half time step implicit in the x -direction, equations (4) and (5) are 
given as 


uP+1/ 2 At 
U i,3 , At 

u n +l . + 

[i , At 1 

uP+i . 

uft 1/2 At 

At 

4 Ax 2R(Ax) 2 

i-l,3 

R(Ax) 2 

1,3 + 

4 Ax 

2R(Ax) 2 


n+i 

U i+l,j 


At 
4 Ax 




) + 

vf-tl/ 2 At 
l >3 ( At 

u u-( 2 + 

\ At 

) + 

4 Ay 2R(Ay) 2 

R(Ay) 2 


uP+1/ 2 

1,3 


At 


vP^l/ 2 At 
i,3 


2R(Ay) 2 4 Ay 


u?tl ( 2 

i,3+l 


( 10 ) 
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and 


u^ 1 * 4 / 2 At 
1,1 , At 

n+1 

[i, .At ; 

n+i 

u P"t^/ 2 At 
1,1 At 

4 Ax 2R(Ax) 2 

Vl,j + 

R(Ax) 2 

V i,j + 

4 Ax 2R(Ax) 2 


.n+1 

1+1,3 


At Ln 


4 Ay 


K i+ i-p"i-i) + 



v n+l/2 ££ 



\ + 

l >] , At 

n+l/2 

It At 1 

j + 

4 A y 2R(Ay) 2 

v i,j-l + 

R(Ax) 2 


n+l/2 
1 ,] 


At 


yPt 1 / 2 At 
v i,3 


2R(Ay) 


2 4 Ay 


n+l/2 

i>j+l 


(ID 


Note that a simple linearization of the convection term is used by evaluating the velocity 
at the previous time step. Since there is interest in only the steady -state solution, this 
assumption should not have any effect bn the results presented herein. 

Equation (7) written implicitly in the x -direction and then in the y-direction is 


— — — pH+,1 • + 

(Ax )2 P i-1,] 


1 + 


At 

(Ax) 2 


f;k+l At -k+1 -Af 1 (n n+1 iiH+1 Vi ' 1 fv I1+ l , 
P U (Ax) 2 i+1 ’i [rat H+i,] u i-i,jJ + 2A^lVi+i 




At L 


(Ay)' 


(pu-i + pf, 


U+l 



( 12 ) 


and 


At -pKt 2 1 + 
2 l d-l 


(Ay) 


fl + At' 

6 k+ 2 

- At p^t 2 - At 

(Ay) 2 

P U 

(iy) 2 P ‘ ,J+1 ‘l 


1 /(.n+1 „ n+1 ^ . 1 /„n+l 

Ax\ i+1 ’j " 2 Ay\ i »i +1 


- v- 1 ,) 
- 1 / 


At /+k+l , -k+1 

+ ^ Pi+1 ’ j+Pi “ 1 > j 


2pM + pH 1 
p i,J / 


(13) 
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where k represents the iteration and the u and v velocities from the last time step 
are used. The boundary conditions for equations (8) to (11) are 

u(i,l) = u(l,j) =.u(i max ,j) = 0 u (Umax) = 1 


v(i,l) = v(l,j) = v(i max ,j) = v(i,j max ) = 0 


The pressure correction p was computed at the boundaries by three methods. The 
first was to solve equations (12) and (13) on the boundaries. For this computation, 
reflection of the node points was used. The second method was to use three-point 
extrapolation utilizing the information that the normal derivative of p is 0. This 
extrapolation is 

Pi=4( 4 Pi + l -Pi+2) (14) 

The third method is to simply set the first interior pressure equal to the boundary value. 
Whether these assumptions are first- or second-order accurate is not material since p 
is zero at steady state. 


APPLICATION 


The procedure for solving equations (8) to (13) is as follows: 

(1) With the initial conditions, u and v are computed first by equations (8) 
and (9) and then by equations (10) and (11). With the velocity of the upper moving sur- 
face equal to 1, the greater change is in the vertical direction and thus the reason for 
the order of solution. 

(2) These velocities from equations (10) and (11) are used in equations (12) and (13) 
to compute p. These two equations (eqs. (12) and (13)) are iterated until the current p 
at each node is within 10 of the corresponding previous value. 

(3) The velocities u and v are then updated to correspond to the new p by 


u Ptl = U P+ X 

u i,J U LJ 


,n+l _ v n+l 

y v u 


+ 2 Ax(pi+l,j 

+ 2~Ay(py+l 
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(4) These values of u and v at time level n + 1 are compared with u and v 
at time level n. If values are within 10 of one another, convergence is accepted. K 
not, p is combined with the previous pressure p for a new value. The pressure cor- 
rection p is subtracted from p at each new time step. .Adding p gives incorrect 
trends in the results. 

(5) The velocities u and v are recomputed in the momentum equations and the 
procedure is continued until convergence is obtained. 

(6) Before recomputing p the initial value of p is always reset to zero. At 
convergence, p should be near zero, so that continuity is satisfied. 

RESULTS AND DISCUSSION 

Solutions are presented for three test cases. For case 1, R = 100, At = 0.01, 
and p was calculated at boundaries with equations (12) and (13). Case 2 is the same 
except R = 10. For case 3, R = 100^ At = 0.005, and p was calculated at boundaries 
with equation (14). 

Figures 2(a) and 2(b) are plots of u and v, respectively, through the vortex cen- 
ter for the previously mentioned cases. The location of the vortex center for case 1 is 
very close to that shown by Mills (ref. 1). For case 3, the vortex center is shown to be 
slightly higher but at the same axial location compared with the results of Mills. This 
is possibly due to slightly different pressure gradients in the two calculations. For 
case 2, the vortex center is at about the same height as for case 1, but the axial location 
has been displaced to the left and the profile is more symmetric. Also the slope of the 
u -profile is lower at the upper boundary. These results can be attributed to greater dif- 
fusion at R = 10. 

In figures 3(a) and 3(b), the pressure coefficient for case 1 through the vortex cen- 
ter is presented in the normal and axial direction, respectively. From the wall (y = 0) 
to y = 0.5, there is good agreement between the present results and those of Mills. 

(See fig. 3(a).) From y = 0.5 to the moving boundary (y = 1.0), there is poor agree- 
ment. Consequently, the agreement between pressure coefficients of the present work 
and those of Mills through the vortex center (y = 0.75) in the axial direction is likewise 
poor. The reason for this disagreement is not known. 

Figures 4(a) and 4(b) are velocity vector diagrams for case 1 and case 2, respec- 
tively. The vectors are of equal length since a characteristic length was used as the 
reference. This was done to amplify the velocities at the cavity bottom. For case 2 
(fig. 4(b)), the possible formation of vortices is indicated in the bottom left and right 
corners. This result is shown in reference 1. When either of the results of figure 4 
are superimposed on the results of Mills, a very good comparison is obtained. 
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A note may be made of the pressure gradients through the cavity. In figure 3, the 
pressures are generally in poor agreement with published results. However, consider- 
ing the good agreement obtained for the velocities in figures 2 and 4, it is assumed that 
the pressure gradients computed in the present analysis must be reasonably well 
predicted. 

Other conditions were considered in investigating this problem. Solving for p on 
the boundaries by the third method discussed previously (e.g., p(l,j) =■ p(2,j)) resulted in 
values for u, v, and p which were not in good agreement with results obtained by 
computing p at the boundaries by equations (12) and (13) or by equation (14). A con- 
vergence criterion on u and v of 10 "4 rather than 10 gave much better agreement 
with results of reference 1. Convergence could not be obtained if accuracy on p was 
raised from 10-3 to 10-4. 

For any of the applications, the values of u and v are not iterated to a required 
convergence at each time step. Convergence is only tested for the values at the current 
and previous time steps. This type of time marching would only produce convergence at 
steady state, and the transient results are not accurate. 

Since most of the documented results on this flow problem are presented in terms 
of stream -function — vorticity variables, an overall comparison of the present results 
computed in primitive variables is difficult. However, where at least a cursory com- 
parison is possible, slightly lower levels for u and v are computed by the method 
discussed in this paper, and a poor comparison for the pressures is indicated. This 
implies that the results obtained by the simplified approach presented herein are not so 
accurate as the results obtained by a more rigorous analysis. 

CONCLUDING REMARKS 

An alternating -direction implicit (ADI) finite -difference procedure has been 
employed to compute the two-dimensional viscous incompressible flow generated within 
a rectangular cavity by a moving surface at the cavity open end. The pressure through 
the cavity was computed by the Spalding SIMPLE procedure. 

Although the present method does not predict axial and normal velocities as high 
as indicated by the more exact methods, the values of the velocities through the vortex 
center are reasonably predicted and the overall velocity vector diagrams appear satis- 
factory. The pressure levels are generally in poor agreement with other predictions. 
However, the pressure gradients must be reasonably well predicted because of the accu- 
rate results obtained for the velocities. 


Ill 
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(b) Profile in x -direction. 
Figure 3.- Concluded. 
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Figure 4.- Velocity vector diagrams. 






9. CUBIC SPLINE SOLUTION FOR THE DRIVEN CAVITY 

Stanley G. Rubin* and Randolph A. Graves, Jr. 

Langley Research Center 

SUMMARY 

A cubic spline procedure for solution of the vorticity— stream -function system for 
flow in a driven cavity has been described. The spline results were compared with finite 
difference results and were found to be more accurate, even near the boundaries. A 
nonuniform grid was also used with increased mesh density in regions where diffusion 
was important. Considerable advantages of the spline procedure with a nonuniform grid 
were indicated in regions of large gradients. 

INTRODUCTION 

The present paper describes a cubic spline procedure for the solution of the driven 
cavity problem. (See paper no. 1 by Rubin and Harris for details of this problem.) A 
vorticity— stream -function system is considered. A finite -difference discretization is 
used for the marching or timelike direction. Unlike a finite -element or Galerkin proce- 
dure, there are no quadratures to evaluate and the coefficient matrix is tridiagonal. The 
spline approximation is second-order accurate even with relatively large variations in 
the mesh. Moreover, for the convective terms (first derivatives) the spline procedure is 
third-order accurate with a nonuniform mesh and fourth-order accurate with a uniform 
mesh. A detailed discussion of splines and spline procedures for fluid dynamics can be 
found in reference 1. 

SYMBOLS 


a,b end points of interval 

hi,kj mesh widths; Xj - X[_i and yj -yj_i, respectively 

L^, Mi spline second derivative in y- and x-direction,, respectively 
#i,mi spline first derivative in y- and x-direction, respectively 


^Visiting professor at Old Dominion University, Norfolk, Va. Now professor of 
aerospace engineering, Polytechnic Institute of New York, Farmingdale, N.Y. 


119 




N number of nodes in interval [a,b] excluding boundaries 

R Reynolds number 

cubic spline function 

t time 

At time step 

u,v axial and normal velocity component, respectively 

X[ nodal points 

x,y coordinate in axial and normal direction, respectively 

£ vorticity 

9 index of finite -difference scheme; 0 for explicit, 1/2 for Crank-Nicolson, 

1 for implicit 

r fictitious time 

At fictitious time step 

if/ stream function 

Superscripts: 

n index of time step At 

s index of fictitious time step At 

Subscripts: 

i,j index of nodal points in x- and y-direction, respectively 

max maximum 
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t differentiation with respect to time 

x,y differentiation with respect to x or y 

SPLINE FORMULATION 
Basic Spline Theory 

Consider a mesh with nodal points (knots) xj such that 
a = x o < x i < x 2 • • • < x N < x N+1 = b 

and with 


h^ = x^ - Xi_i > 0 

Consider a function u(x) such that at the mesh points Xi, 

u(xi) = Ui 


The cubic spline is a function S^(up£) = S^(x) which is continuous, together with its 
first and second derivatives, on the interval [a,b], corresponds to a cubic polynomial in 
each subinterval Xi_i lx S Xi, and satisfies S^uipq) = Ui. 

If S^(x) is cubic on [xi.^xj, then in general 


S’A(x) = Mi. 



where 


Mi = S^(xi) 


Therefore, 


S A (x) = Mi.x 


<X| ‘ X>3 I M- <X ' X '- 1>3 l 1 

L . M i-i h i 2 ^ 

IX; - X j 

( Mihi 2 ^ 

1* - x i-l 

6hi Ml 6 hi 

f 1 " 1 ‘ 6 ) 

^ +l 

V Ul 6 j 

' hi 


(1) 
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The unknown derivatives are related by enforcing the continuity condition on 
S^(x). With S’(Xj-) = mj- on [x^^xj and xj = m* on [x if x i+1 ], 

mf = mf = mj 
For i = 1, . . N, 


hi hi + hiil llij.1 


M; 


i+1 


Uj+1 - u i _ Uj - uj-l 
hi+1 " hi 


( 2 ) 


h i .. hi ^ - u 

m i = y M i + T Mi - 1 + ^ 


i-1 


(3) 



Mi 


hj+1 

6 


M i+ i + 


u i+l ■- u i 
hi+1 


Equation (2) lead to a system of N equations for the N + 2 unknown Mi. 


(4) 


Splines for Solving Partial Differential Equations 

If the values ^ are not prescribed but represent the solution of a quasi -linear 
second-order partial differential equation, such as 


u t = f ( u > u x» u xx) 


then an approximate solution for Ui can be obtained by considering the solution of 

(u t )i = f(ui,mi,Mi) 


where the time derivative is discretized in the usual finite -difference fashion; that is, 
n+1 n 

^ = (i . e)t n + ef n+1 

At 

where 9 = 0, 1, and 1/2 for explicit, implicit, and Crank -Nicolson finite -difference 
schemes, respectively. The physical time t equals n At (At is the time increment 
at each step). 
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For equations with two space dimensions such that 


U t = f^Ux.Uy.Uxx.Uyy) 

a spline -alternating -direction -implicit (SADI) formulation is developed. The two-step 
procedure, with quasi -linearization or some other iterative process used for nonlinear 
terms, is of the following form: 

Step 1 


P+i/2 = 




(5a) 


Step 2 




(5b) 


where and Ly are the spline approximations to (u y ). . 

respectively. ^ 


and 


( u yy)i,j 


INCOMPRESSIBLE FLOW IN A CAVITY 


Equations 

For the driven cavity the governing equations, in terms of vorticity and stream 
function, are 


^xx + ^ yy ^ 


^t + u ^x + v ^y r(^xx + *>yy) 




( 6 ) 


where 4 * is the stream function, £ is the vorticity, and u = i^y and v = -i^ x are 
the velocities in the x- and y-directions, respectively. For all calculations the initial 
conditions are ^(x,y) = 0 and C(x,y) = 0. 

Solutions are obtained by the iterative SADI procedure. The SADI system repre- 
senting equations (6) is given in two steps for both stream function and vorticity. 
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Stream function. - For step 1, 


.n+l, s+ 1/2 .n+1,8 at 

^i,j 2 


hT^-kT 8 -* c 


n+l 

j 


(7a) 


For step 2, 


.ai+l,s+l _ wji+1,8+1/2 


♦VHT ”*' 2 


+ 



n+l, s+l 


»>n+l 

s i,j I 


(7b) 


where At is a fictitious time step and t = s At. Solutions for equations (6) are 
obtained as the steady-state limit (t — <*) of equations (7). The notation Ly and My 
denote the spline approximations to a^A/ay^ and a^A/ax^, respectively. (The super- 
script ijs implies that A = \p.) First derivatives i//y = u and \f/ x = -v are repre- 
sented by and m'y, respectively. 


Vorticity. - For step 1, 



For step 2, 




n+l 


* 


. n+1/2 


+ R- 


n+1/2 


Hi) ' -- 1 Ki 


n+l 


(8b) 


The bar over superscript <// in the spline derivatives denotes an average of values at 
time steps n and n+l. The £ superscript denotes C spline derivatives. 

Numerical Solution Procedure 
The iterative procedure is as follows: 

(1) With i/'y and £y given either as initial conditions or at time n At, all the 
i// spline derivatives are determined from equations (1) to (4). On the vertical surfaces, 
on the horizontal boundaries, Lj’H = .. 


124 



n I 1 

(2) The vorticity j is obtained with the SADI technique outlined by equa- 
tions (5). At the boundaries, is found from an expression similar to equation (3) 

or (4). At the upper moving wall (w), 


„<// t _ k i,w L ^w k i,w L i^w-l ^i,w " ^i,w-l 
n,w ' 1 3 + g + jT 

O O K 1,W 


where kj = yj - yj _ x . With ^ i>w = 0 and L^ w = 


w> 




3 L i,w -1 3 4 * ^i,w-l 
‘ + , 2 


w 


k i,w ^ 


n,w 


Similar relations can be derived for the three stationary walls. Boundary values for 

. and L? . are obtained from equations ( 8 ) evaluated at the surface. For the mov- 
i,j bJ jo I/O 

ing wall, equations (3) and (4) are used to eliminate m? j. In addition, is eva l- 

uated with the three -point formula; that is, 


? „ +w i/ 2 . 3 ^w^. 6 |,w ;ji,w + 0 [ (4t)2 ] 


(3) The vorticity is used in equations (7) and the SADI procedure is applied 

1 >J 

over the fictitious time r until a converged solution, to any specified tolerance, is 
obtained. 

(4) If only the steady -state solution is required, the calculation proceeds to the 
next time step n + 2 by returning to step (2) with n replacing n + 1. The spline 
derivatives for \}y have already been determined in step (3). 


If an accurate transient solution is required, the calculation proceeds to step (2) 
with i/'j j and all spline derivatives of »//' replaced by averages over the n and n + 1 
time steps. Then CPt* is recalculated. This process continues until convergence is 
reached. 


bJ 


Results 

Calculations are presented for a square cavity with R = 100. Comparisons are 
given with finite -difference calculations in both divergence and nondivergence form as 
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obtained from the results of Morris (paper no. 5) and Smith and Kidd (paper no. 6) in this 
volume and of Bozeman and Dalton (ref. 2). 

The results are presented in tables I to IV and in figure 1. In all instances, the 

values of t// max and the vorticity ? wa ll at the midpoint of the moving wall are 

depicted. For these quantities, comparisons between the spline and finite -difference 

solutions are possible even when the grid alinements differ. In addition, the distributions 

of I// and £ for the spline solutions and, in several cases, for the finite -difference 

solutions are presented. Figure 1 depicts the horizontal velocity component uy or 

i'P . along a vertical line passing through the vortex center. 
bJ 

The values of >/ / max and are shown on table I for a variety of grids. 

Also included is a limiting solution obtained by Richardson extrapolation from the two or 
three calculated values of each procedure. It is evident that the divergence finite - 
difference solution is more accurate than the nondivergence result; however, the spline 
solution, which is obtained in nondivergence form, appears to be even more accurate than 
the divergence finite -difference result. For example, the value of i^ max obtained from 
the spline calculation with a uniform grid of 15 points in each direction is about 1 percent 
higher than the extrapolated value; the 17 -point divergence finite -difference result is 
about 4 percent lower than its extrapolated value. The nondivergence 15-point finite- 
difference result is low by about 12 percent. These results seem to reflect the higher 
order accuracy of convection terms in the spline procedure. In the vortex core region 
the flow is inviscid dominated. However, near the moving wall where diffusion is impor- 
tant, the vorticity results appear to show a similar trend; the spline values are always 
somewhat more accurate than the divergence finite -difference solutions. 

Another interesting result is shown on table 1(c). The velocity at the first grid 
point away from the upper left boundary is depicted. With 15 points the spline result has 
a sign opposite that obtained with the nondivergence finite -difference method. With a 
finer grid the finite -difference solution changes sign, so that once again the spline pro- 
cedure prevails. Unfortunately, the more accurate divergence finite-difference solutions 
were obtained with a slightly different grid, so that a direct comparison is not possible. 
However, a change in sign with mesh reduction is observed for the velocity in the corner 
region. An interpolation procedure is used to estimate the value at the desired location. 
These results are also given on table 1(c). The extrapolated limit closely approximates 
the solution obtained with splines. The velocity profiles through the vortex center are " 
shown in figure 1. The solutions for i// and £ are given in table II. 

A spline solution was also obtained with a 19-point nonuniform mesh. In the cen- 
tral region of the cavity, hj = 1/14, as with the 15-point mesh; however, near the bound- 
aries there is some grid realinment to increase the mesh density in the surface boundary 
layers. The increased accuracy near the boundaries, where diffusion is most important, 


126 



leads to a solution that is almost as accurate throughout the entire flow domain as the 
29 -point results. The improved accuracy of this 19-point solution is seen on table I 
where i// max and £ wa u are indicated. These results imply the considerable advan- 
tages of the spline procedure with a nonuniform grid in regions of large gradients. With 
the spline procedure the accuracy of the second -order diffusion terms is enhanced in 
domains where diffusion effects are significant. In inviscid regions the fourth-order 
accurate convection terms are dominant and mesh reduction is not as important. The 
improved resolution of the corner vortices is seen in the \p and £ distributions of 
table HI; the comparisons with the 65-point divergence finite-difference solutions 
(table IV) are reasonably good. 


CONCLUDING REMARKS 

A comparison between results of the spline -alternating -direction -implicit (SADI) 
and finite-difference procedures indicated that the spline solution for flow in a driven 
cavity was more accurate, even near the boundaries. A nonuniform grid was also used 
with increased mesh density in the boundary layers where diffustion was most important. 
The spline solution with a 19 -point nonuniform grid was almost as accurate throughout 
the entire flow domain as the spline solution with a 29-point uniform grid. Thus, consid- 
erable advantages of the spline procedure with a nonuniform grid are indicated in regions 
of large gradients. 
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TABLE I. - RESULTS FOR SQUARE CAVITY WITH R = 100 


(a) Vorticity at center of moving wall 


Calculation method 

Points 

Vorticity at center 
of moving wall 

Spline 

15 x 15 

7.1376 

Spline 

29 x 29 

6.6876 

Extrapolated spline 


" 6.5376 

Spline (unequal spacing) 

19 x 19 

6.2970 

Finite difference 

15 x 15 

8.9160 

Finite difference 

57 x 57 

6.6960 

Extrapolated finite difference 


6.5480 

Finite difference 3 

17 X 17 

7.3755 

Finite difference 3 

33 X 33 

6.7653 

Finite difference 3 

65 x 65 

6.6091 

Extrapolated finite difference 3 


6.5567 


a Divergence form. 


(b) Maximum stream function 


Calculation method 

Points 

Maximum 
stream function 

Spline 

15 x 15 

-0.10529 

Spline 

29 x 29 

-.10432 

Extrapolated spline 


-.10399 

Spline (unequal spacing) 

19 x 19 

-:10472 

Finite difference 

15 x 15 

-.08742 

Finite difference 

57 x 57 

-.10128 

Extrapolated finite difference 


-.10220 

Finite difference 3 

17 x 17 

-.09867 

Finite difference 3 

33 x 33 

-.10213 

Finite difference 3 

65 x 65 

-.10318 

Extrapolated finite difference 3 


-.10355 

Reference 2 

51 x 51 

-.10316 


divergence form. 


(c) Corner point velocity u 


Calculation method 

Points 

Velocity u 
at corner 
(a) 

Spline 

15 x 15 

-0.13230 

Spline 

29 x 29 

-.10036 

Extrapolated spline 


-.08971 

Finite difference 

15 x 15 

.05730 

Finite difference 

57 x 57 

-.06615 

Extrapolated finite difference 


-.07438 

Interpolated finite difference 15 

17 x 17 

.02079 

Interpolated finite difference 15 

33 x 33 

-.05189 

Interpolated finite difference 15 

65 x 65 

-.07 560 

Extrapolated finite difference 15 


-.08399 


a Locations: 0.07143 down from moving surface; 0.07143 in from 
left surface. 


^Divergence form. 






TABLE n.- VORTICITY AND STREAM FUNCTION CALCULATED WITH SADI PROCEDURE WITH EQUALLY SPACED 15 X 15 GRID 
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TABLE IV.- VORTICITY AND STREAM FUNCTION CALCULATED WITH DIVERGENCE FINITE -DIFFERENCE FORM WITH EQUALLY SPACED 65 x 65 GRID 


(a) Stream function 
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Figure 1.- Comparison of calculated u through pc 



Nondivergence Form 
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Finite Difference, 57 x 57 points 
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Finite Difference, 65 x 65 points 


J I I I I 1 , I 

4 .5 .6 .7 .8 .9 1.0 


of maximum for R = 100. 
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